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I. INTRODUCTION 


Noise generated by the turbulent mixing of jet engine exhaust flow with ambient air has 
become a problem of greater concern with the advent of plans for a high speed civil transport, HSCT. 
The viability of the HSCT depends on being able to satisfy airport noise regulations. Efforts to 
develop noise suppression techniques solely on the basis of experimental work are extremely expensive 
and can easily miss directions that could result in significant advances. Computational approaches 
to date rely heavily on empirical results and/or greatly simplified approximations of the governing 
equations. Because much of the basic flow and acoustics is not well understood, success in computing 
noise for one nozzle and set of flow conditions does not guarantee success for other geometries which 
may employ different physical concepts for suppression. 

The work reported here is based on a fundamental study of both the time-dependent fluid 
mechanics and acoustics of turbulent jets in three dimensions. The emphasis will be on developing 
procedures for the computation of far field noise of interest for FAA certification and community 
noise studies. The methods can also be used for near field sound problems of interest, for example, 
for airframe fatigue evaluations. 

The time-dependent approach is essential to understanding the basic mechanisms. The time- 
dependent Navier-Stokes equations are solved using a spectral element technique in conjunction with 
renormalization group theory, RNG, to treat the small scale turbulence. RNG facilitates computation 
time reductions by orders of magnitude since it properly models the effect of the small scale 
fluctuations on the large scale fluctuations without having to explicitly compute their behavior. This 
can be applied for the acoustics problem because it is the larger scale fluctuations that are the most 
efficient for sound production. While the jet flows studied here are for the simplest possible nozzle 
geometries the numerical methods incorporated are perfectly general, and could thus be applied to 
much more complex nozzle shapes. One of the major benefits of the spectral element method is its 
ability to treat complex geometries as would be encountered in noise suppressors. 

Figure 1 outlines some methods that can be considered for obtaining the far field noise 
spectrum. If the fully compressible time-dependent Navier-Stokes equations are solved, then the 
compressible pressure field will be automatically computed in the ambient region surrounding the 
turbulent flow. This pressure could then be extrapolated to any desired far field location using 
standard techniques based on the solution of the acoustic wave equation. Compressible flow 
computation methods will be discussed, but acoustic results will be generated with incompressible 
flow information in this report. 

The second approach, utilized here, is to solve the incompressible Navier-Stokes equations and 
use the time-dependent turbulent flow fluctuations as volume source terms to force the wave 
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equation. Acoustic solutions are obtained by a finite difference method. Thus, compressibility enters 
the problem only through the acoustic wave equation. Far field solutions can be obtained by 
extrapolation of the acoustic field at a distance of a quarter wavelength from the jet flow 1 as 
described in Appendix A. 

This report begins with a description of the flow field computation methods and the general 
acoustic analogy approach to finding the sound field in terms of the Lighthill, Lilley, and Phillips 
equations. The finite difference solution to the wave equation is obtained, tested, and generalized 
to include the effects of convection by a time-dependent flow field. Methods are then presented for 
extrapolating near field pressure to far field locations. Flow results are presented for a two- 
dimensional jet, and both flow and noise results for a three-dimensional jet. Finally, 
recommendations are made for improving the analysis procedures. 

We have attempted to discuss both theoretical and numerical aspects of the problem with 
enough detail so that various aspects of the analysis could be utilized by the reader in other situations. 
A computer operations guide is included as Appendix C of this report and should aid in allowing one 
to change the various conditions and parameters found in the combined flow and acoustics programs. 
A separate manual for running the flow computation pre- and post-processors on workstations is 
being provided to the contract monitor. All of the software that we have used to obtain flow and 
acoustic results has been either transferred to an account accessible to NASA personnel on NASA 
computers or delivered to NASA. Restrictions on the use of the code are discussed in Appendix C. 

We greatly appreciate the advice of Prof. Juan Ramos of Carnegie Mellon University who 
consulted with us during the early phases of our work and the advice and encouragement of members 
of the Naval Research Laboratory’s Computational Physics Group, especially Drs. Elaine Oran and 
Fernando Grinstein. 


II. FLOW COMPUTATIONS 


The assumption made in this paper, based on other work, is that the turbulent fluid mechanics 
is well represented by incompressible computations over a significant range of Mach numbers. 
However, Mach number effects are extremely important in the acoustic problem since the self- 
cancellation tendency of multipole sources is reduced both by convective motion and a shift to higher 
frequency fluctuations as jet velocity increases. 

The computational approach chosen uses the RNG turbulence model and a Navier-Stokes 
equation solution routine called NEKTON. From the practical point of view both of these features 
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complement each other and provide a solution method that can obtain the information needed for 
noise computations at high Reynolds numbers within reasonable computation times for studying the 
basic mechanisms. 


A. RENORMAI.T7.ATTON GROUP METHODS 

RNG methods involve systematic approximations to the full Navier-Stokes equations that are 
obtained by using perturbation theory to eliminate or decimate infinitesimal bands of small scale 
modes, iterating the perturbation procedure to eliminate finite bands of modes by constructing 
recursion relations for the renormalized transport coefficients, and evaluating the parameters at a 
fixed point in the lowest order of a dimensional expansion around a certain critical dimension. The 
decimation procedure, when applied successively to the entire wavenumber spectrum, leads to the 
RNG equivalent of full closure of the Reynolds averaged Navier-Stokes equations. The resulting 
RNG transport coefficients are differential in character as opposed to ad hoc algebraic coefficients 
of conventional closure methods. All constants and functions appearing in the RNG closures are fully 
determined by the RNG analysis. 

In essence, the RNG method provides an analytical procedure for eliminating small scales 
from the Navier-Stokes equations, thus leading to a dynamically consistent description of the 
large-scales. The formal process of successive elimination of small scales together with re-scaling of 
the resulting equations results in a calculus for the derivation of transport approximations in turbulent 
flows. The governing equations for the flow motion are the Navier-Stokes equations: 


3v- 3v- 

+ V: — 

dt 3 3xj 


1 dp + d_ 
p 3Xj dXj 



along with the incompressibility constraint 


5v i 

3 Xi 


0 


( 1 ) 


( 2 ) 


Here v is the total viscosity defined as v = + i/ x (the sum of the molecular and turbulent viscosity, 

respectively). 
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For large eddy simulations the eddy viscosity v is obtained from the following relation: 


1 + H 

a£ A 1 * * 4 - cl 



(2x) 4 i/J 

L -d 



( 3 ) 


where H(x) is a ramp function defined by H(x) = x for x > 0 and H(x) = 0 otherwise, and A is the 
width of a suitably chosen Gaussian filter equal to twice the computational mesh size. The mesh size 
must be small enough to resolve scales that lie in the inertial subrange of the turbulence spectrum. 
The constants a = 0.120 and C = 100 are derived in Ref. 2. The mean dissipation rate e can be 
expressed entirely through the resolvable field as 


e 


v 

2 



3Xj 3X; 


^ 2 


( 4 ) 


The total viscosity is then obtained by solving an algebraic cubic equation at every node of the 
computational domain at each time step in terms of the dissipation at the previous time step. 


B. NF.KTON CODE 

NEKTON solves the steady and unsteady three-dimensional primitive variable Navier-Stokes 
equations using preconditioned iterative (multigrid) solvers and global finite-element spatial 
discretizations to achieve high accuracy.® The code is being extended to include the effects of 
compressibility. The NEKTON code has advantages over other computational approaches not only 
in efficiency, but also in generality, robustness, and its potential for parallel computation. 


1. Time Stepping and Compressibility 

The NEKTON code uses a combination of high-order explicit and semi-implicit 

unsteady time stepping procedures that capture the flow properties without the penalty of (diffusively 

restricted or sound speed restricted) small time steps. The nonlinear terms of the Navier-Stokes 
equations are handled with third-order accuracy in time in order to achieve both high relative 
stability and accuracy. Together with the high spatial accuracy of the global finite— element 
formulation, the NEKTON code minimizes numerical dissipation and diffusion errors. 
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NEKTON works efficiently across a broad range of Mach numbers up to about M = 
2 for one-dimensional time-dependent flows. For higher-dimension flows it is restricted to slow time 
variations and is thus not suitable for studying turbulence. 

A three-dimensional, compressible, time-dependent code has now been written to treat 
turbulent jet flows. This code, suitable for research studies, has been delivered to NASA without any 
restrictions on its usage since it is not a part of the general NEKTON code. It is not yet ready for 
production runs. Restrictions on the use of the NEKTON code are discussed in Appendix C. This 
and other aspects of compressibility are discussed separately in Appendix B since no acoustic analysis 
was performed in conjunction with the compressible code. 

The current NEKTON does not suffer from severe time step restrictions, even if the 
Mach number is relatively small. At small Mach numbers an implicit procedure for density evolution 
at each time step is based on a pressure-velocity decoupling and inner-outer Richardson iteration that 
can be proved to result in error decrease by a factor of ten at each iteration. This result should be 
compared with older SIMPLE 4 iterative procedures that can require thousands of iterations, and are 
very problem-dependent in terms of convergence history. The detailed formulation in NEKTON 
yields minimal boundary errors associated with the pressure-velocity decoupling and yet enables 
efficient solution of the Navier-Stokes equations in the complex geometries of engineering interest. 

2. Spatial Accuracy. Discretization Error, and Fine Grid Resolution 

The NEKTON code uses global finite-element spatial discretizations, in which the 
computational domain is first broken up into macro brick-shaped global finite-element (spectral 
element) subdomains. (The individual macro-bricks will have curved sides to follow solid boundaries, 
but are topologically parallelepipeds.) The dependent and independent variables in each element are 
then written as high-order polynomial expansions in terms of tensor-product Legendre-Lagrangian 
interpolants, and variational projection operators are used to generate the global discrete equations. 
The global finite-element method is applicable in virtually arbitrary geometries due to the 
macro-element decomposition. NEKTON may be run to solve the unsteady full Navier-Stokes 
equations or the steady or unsteady linearized (low Reynolds number) Stokes equations, with or 
without heat transfer effects. 

The accuracy of the global finite-element method may be summarized as follows: The 
convergence rate of the method depends only on the smoothness of the underlying solution, with 
infinite-order convergence obtained for analytic solutions. In contrast to finite-difference or h-type 
finite-element solutions in which doubling the spatial resolution decreases the error by a fixed factor 
(4 for a second-order scheme), the error in NEKTON decreases exponentially fast. Furthermore, a 
wide variety of studies have shown that to achieve engineering accuracy of several percent the global 
finite-element method requires at least an order-of-magnitude less resolution in two-dimensional 
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flows than do difference or h-type element methods. Once engineering accuracy has been achieved, 
further increase of the spatial resolution in NEKTON leads to exponentially small errors because of 
the high-order character of the method. The global finite-element discretization of the convection 
operator results in minimal numerical diffusion and dispersion, a fact which is of critical importance 
in the solution of difficult laminar or turbulent flows. Also, the global finite-element method offers 
good resolution of boundary layers due to both the intrinsic clustering of high-order polynomial 
interpolation and the flexibility of an elemental decomposition. The accuracy of the global 
finite-element discretization is relatively insensitive to skew or disparate-length-scale elements due 
to the fact that the discretization is variational in nature. Also, in contrast to difference or h-type 
element methods, errors due to geometrical patching or grid-resolution changes are small with the 
method capable of handling rapid changes in resolution across interfaces. 

NEKTON’S efficiency is further enhanced by the use of multigrid accelerated iterative 
elliptic solvers (for the pressure and viscous terms in the Navier-Stokes equations) which are superior 
to direct solvers as regards both storage and operation count. Furthermore, these elliptic solvers are 
well-suited to medium-grain parallelism, and have already been implemented with near-unity parallel 
efficiency on the Intel iPSC Hypercube, a prospect that may be of greater interest as these parallel 
machines develop further in later years. 

The work reported here will be limited to simple geometries: a two-dimensional slot 
nozzle with two-dimensional turbulence and a round nozzle with three-dimensional turbulence. The 
two-dimensional case may not be of direct interest for jet noise, but the inclusion of a specific nozzle 
lip geometry and its effect on the development of the turbulent field is. 

We considered the direct numerical simulation of a two-dimensional jet at a Reynolds 
number of 300 based on the jet exit dimension and the maximum velocity. Thus, no modeling of the 
turbulence is employed. The jet exit nozzle is represented by two parallel semi-infinite plates, or 
fins, that protrude one unit length, L, into the computational domain, are a distance of L apart, and 
are each 0.2 L thick. The jet velocity profile in the nozzle is parabolic with a maximum value of 
unity. A wall boundary condition, i.e., zero velocity, is imposed on the upstream boundary of the 
rectangular computational domain above and below the fins, periodic conditions are specified on the 
sidewalls located at ± 22 L from the jet centerline, and an outflow condition is specified at the 
downstream wall, a distance of 64 L from the upstream boundary. A highly viscous damping region 
exists just upstream of the outflow boundary to reduce the effect of reflections caused by the 
interaction of the large scale eddies with that boundary. 

An instantaneous vector velocity plot of the flow field is shown in Fig. 2a at Re = 300. 
A similar plot, Fig. 2b, was obtained utilizing an RNG large eddy simulation model at Re = 10,000 
for the same equal fin case. The major difference between the high and low Reynolds number cases 
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are the much larger fluctuating velocities at high Reynolds number within 5 to 10 jet diameters of 
the nozzle exit. At much greater distances the results are qualitatively similar. 

The rest of this report will focus on the round nozzle simulation. The gridding pattern 
for this is shown in Fig. 3a in a plane normal to the jet axis with seven nodal points per element. 
Figure 3b shows the three-dimensional element boundaries. Other details of the grid and flow will 
be discussed later. 


III. ACOUSTIC COMPUTATIONS 

A. NOISE SOURCE THEORIES 

The most well-known method is that due to Lighthill, 5 which results in an integral equation 
solution involving turbulent noise source terms. 

Lighthill’s equation is 


a 2 p _ 2 d 2 p 

3t 2 dxf 

where the turbulent noise source tensor Ty is 

Tjj = P v i v j + 

The last two terms in Ty involving perturbations about the mean in the stress tensor py and the 
density p are considered to cancel each other out by assuming that Py is the hydrostatic pressure and 
thermodynamic fluctuations are isentropic, i.e., dp/d p = c 0 2 . For low Mach number incompressible 
flows we will set p = Pq in the first term of Ty. Thus, it will be fluctuations in the velocity vector V; 
that produce sound. 

Other improvements over the Lighthill equation are found in Phillips’® and Lilley’s^ equations. 
Phillips’ equation treats the problem of convection and refraction of sound. For nearly incompressible 
isothermal flows the source term is nearly identical to Lighthill’s. In Lilley’s equation the explicit 
effect of flow gradients on acoustic propagation is treated in the wave operator. 


d 2 Ty 
dXjdXj 


Pij " c 0^ij 


(5) 


( 6 ) 
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Phillips’ equation is 

D 2 g _ d 
Dt 2 ax i 


da 

<3X: 


au: D 

= 7 + 1 — 

3xj dx { Dt 

+ viscous terms 


a = 
Dt 


_L 

Cp Dt 


In— 

Po 


" + U„ 

dt 


dx u 


( 7 ) 


where S is the entropy and 7 is the ratio of specific heats. 
Lilley’s equation is 


P 3 g 

Dt 3 


D 

Dt 


3U, ^ 

> D 

3xj 3xj Dt 

i D 

1 DS 

7 Dt 

c p Dt 


,2 da 


+ viscous terms 


aU: au k au= 
- 27 — 5 — - — l 
9xj 3x ; 3x k 


( 8 ) 


We believe that either the Phillips or Lilley equations are needed to describe sound production 
in high speed jet flows unless fully compressible methods for solving the Navier-Stokes equations are 
available. Lilley’s equation is very popular for finding the Fourier transform of the acoustic pressure 
fluctuation so that the frequency spectrum can be determined directly. This involves using the time- 
averaged flow field values and often assuming that the flow velocity vectors are parallel and change 
little in magnitude in the mean flow direction. In this case a singularity appears in the differential 
equation and a number of methods of removing that singularity have been proposed. The singularity 
can be thought of as a resonance in which a forcing function excites an oscillator. The simplest 
problem is the linear perturbation of a vortex sheet by a sound wave . 8 Since the Lilley equation uses 
the momentum equation in its derivation, it provides more of a simultaneous solution of the 
momentum and pressure equations as compared to the Phillips equation, which assumes that one has 
an independent means for obtaining velocity. 

The question is which equation is more applicable for our purposes. In a real turbulent jet, 
velocity fluctuations on the order of ± 15% are typical and larger ones are quite likely. For this reason 
one must ask if parallel flow conditions, which are responsible for the resonant singularity, exist in 
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a turbulent flow. If inhomogeneities in the flow field exist, then it may make more physical sense 
to think of the solution in terms of a scattering problem. We note that there is no readily apparent 
singularity or resonance in the time— dependent Lilley equation. For the Phillips equation it is 
important that the velocity fluctuations include the effect of the pressure on them. In our formulation 
the incompressible pressure field is solved for simultaneously with the velocity field in the NEKTON 
code. It would undoubtedly be preferable to utilize a compressible flow code to compute the velocity 
fluctuations. Our viewpoint is that the nature of the turbulent velocity fluctuations does not change 
much from the incompressible case until one reaches a jet Mach number of at least 2.0. The pressure 
ratio needed to achieve such a Mach number is well in excess of any that is anticipated for a HSCT 
engine. From this standpoint both the Phillips and the Lilley equations are equally applicable, but 
we have chosen to use the Phillips equation for practical reasons. The Lilley equation has a more 
complex source term dependent on the products of three velocity derivatives whereas the Phillips 
equation involves the products of two velocity derivatives. Also, the Lilley equation is third order 
in the time derivative which requires another difference equation and quite likely a smaller time step 
increment. 


B. ZERO FLOW CASE 

We next develop time-dependent finite difference solutions of Lighthill’s equation with the 
source term written similar to the source in Phillips’ equation. This formulation would be applicable 
to both the Phillips and Lilley equations for the case of no flow. 

The time derivative in Lighthill’s equation is written in finite difference form as 


p u (x,y,z,t) = (p"j£-2p" jik + 0/(At) 2 


(9) 


and the x coordinate derivative is written as 


Pxx (X,y,z,t) 


( 10 ) 


where the indices i, j, k and n refer to finite difference steps in x, y, z and t, respectively, and similar 
y and z derivatives are formed based on the relations 
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t = nAt 
x = iAx 

y = jAy 

z = kAz 


The finite difference equation is 


+1 


2(1 - h x - h y - h.) P " jik - p-; ik 

+ ^x^i+l,j,k + ^i-l,j,k) 

+ h y^!j*l.k + P U-l,k) 

+ h «(Pu,k.l + ^j,k-i) + s u,k ( At > 2 


( 11 ) 


( 12 ) 


where 


h x = c 2 (At/Ax) 2 
h y = c 2 (At/Ay) 2 
h t = c 2 (At/Az) 2 


(12a) 


and 


5 i,j,k 


dx p 3x q 


p 0 


au p au q 

dx q 3x p 


(13) 


using a strict incompressible representation of the source term. The form of the velocity derivatives 
is the same as in Phillips’ source terms. The quantities h x , h y , h E are the squares of the acoustic 
Courant number C A . 


To date we have been successful in reducing the effect of reflections from the boundaries of 
the computational domain with rather simple acoustic boundary conditions of the form 


— + c 

I <9t dxjsj 

where x N is the component normal to the boundary. An example of the effectiveness of this 
boundary condition is shown in Fig. 4 in which a point source on the centerline of a computational 
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cylinder is turned on for one sinusoidal cycle in time and then turned off. Note that the sound field 
has almost vanished after the wave has struck the absorbing boundaries whereas the wave continues 
to reverberate for the reflecting wall case. 

Examples are shown of some simple test cases with a monopole point source in Fig. 5 and for 
a dipole source composed of two out of phase point sources in Fig. 6. Straightforward analytical 
solutions for both of these cases agree quite well with the numerically computed results. 


C. JF.T FLOW CASE 

For the case of sound wave propagation in a jet flow we have previously attempted to solve 
a convected wave equation, essentially the Phillips equation, using upwind differencing in a fixed 
Eulerian reference frame. Here the partial time derivative in the wave equation (Eq. (5)) without 
flow is replaced by the convective derivative. In the rest of our work here using the Phillips equation 
only the constant speed of sound case will be treated. With a uniform steady flow we encounter 
increasing dissipation of the amplitude of sound waves propagating against the flow with increasing 
Mach number. This is probably due to the low order accuracy of the cross derivative in space and 
time whereas all derivatives are accurate to second order in the standard wave equation without flow. 
For this reason we chose a Lagrangian approach for the flow case so that higher order accuracy could 
be obtained by applying the wave equation operator in reference frames moving with the local flow. 

Methods for solving the time-dependent wave equation in one space dimension for the case 
of a uniform flow were studied with the waves created by a point source varying sinusoidally in time. 
A semi-Lagrangian technique to be discussed next was developed. The term semi-Lagrangian is used 
to describe an approach in which the numerical algorithm utilized is solved in a Lagrangian reference 
frame moving with the local flow velocity viewed in a fixed Eulerian frame. However, after applying 
the algorithm during one time step, one immediately transforms the results back to the fixed frame 
via interpolation and specifies that the coordinates of a new Lagrangian frame for the next time 
interval coincide with the fixed system points. Thus, a new Lagrangian coordinate system is defined 
after each time step so that a given fluid element is followed only during a single time interval. 
Redefining the Lagrangian coordinates after each time step eliminates the need to introduce and 
remove Lagrangian fluid elements and automatically gives results at fixed points where results from 
the NEKTON flow code are obtained. 

In the case when the flow Mach number, M, is zero, the Eulerian and Lagrangian frames are 
identical; a solution is shown in Fig. 7a for the case when there are two acoustic wavelengths in the 
computational domain. 
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When M is nonzero and constant, the two reference frames are related by 

z = z 7 + cMr 


( 15 ) 


t = r 


(16) 


where z and t are the Eulerian space-time coordinates and z' and r are the corresponding Lagrangian 
coordinates. At the current time, grid points i for the Eulerian system coincide with Lagrangian grid 
points j. At the next time step for which the acoustic pressure is to be calculated the j system moves 
a short distance relative to the i points. The solution for the acoustic pressure, p, is obtained in terms 
of values of p at both the current time and the preceding time step. Values at the preceding time step 
are obtained by interpolation of results at that time from the Eulerian to the Lagrangian frame. The 
newly determined Lagrangian results are then interpolated back into the Eulerian frame. Before the 
next pressure value is computed the i and j systems are made to coincide. Thus, the method requires 
continual interpolation back and forth between the two frames. 

Our interpolation scheme uses a Taylor series expansion in space to second order in 
displacement. The first and second derivatives of pressure in the expansion are obtained using finite 
difference representations accurate to second order. The distances between the locations where 
pressures are known and those where they must be interpolated are simply the local speed cM 
multiplied by the time step increment. 

This is shown as follows for a velocity component only in the z direction. The z velocity is 
a function of all three space coordinates and time. The z coordinate is in the fixed Eulerian frame 
and z' is in the moving Lagrangian frame. The expansion is 


p(x,y,z,t) = 


The first derivative on the right hand side of Eq. (17) is evaluated by centered differences and the 
second derivative by the standard expression of Eq. (10) so that the resulting formula is second order 
accurate. The grid spacing Az appears in the finite difference form of the derivatives but the 
displacement, Az', is given by 


p(x,y,z 7 ,t) + (z-z 7 ) — -p(x,y,z / ,t) 

dz' 

*- p(x,y,z / ,t) 


(17) 


dil 


n. 


Az 7 = -(z-z 7 ) = U B At 


(18) 
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Forgetting about the x, y, t dependence, the Eulerian values, given by the left hand side of Eq. (17) 
and denoted next as p E , are written in finite difference form in terms of Lagrangian values p L as 

p E (m) = p L (m)(l - C F )-0.5C F (p L (m + 1) - p L (m - 1)) 

+ 0.5Cp(p L (m + 1) + p L (m - 1)) 


where C F is a Courant number based on flow speed 

C F = Az'/Az = U,At/Az (20) 

Here we don’t want C F to be too large to maintain accuracy in the expansion. Since C F is the product 
of the acoustic Courant number and the flow Mach number, limits on C^, to be discussed later, 
will keep C F in the proper range as long as the Mach number is less than some supersonic value that 
is not too large. 

These interpolation formulas are incorporated into subroutine shftw. The inclusion of only 
the z component for solving the Phillips equation is made from the view that the effect of the flow 
is one of refraction and convection of the sound waves and U E has the largest of the velocity 
component values. If scattering effects are important, then all velocity components should be 
included. This is a straightforward procedure to program, but our decision to not include the other 
components is based on the computer time needed to interpolate the additional flow information. 
New improvements in the computer code will probably allow us to perform interpolations more 
quickly. 

In implementing the above routine we decided to use a local space average value of U E 
obtained as the average of the point in question and the six nearest points surrounding it. This 
averaging is performed in subroutine fmachw which also determines the z component value of the 
Mach number. 


Results are shown in Fig. 7b for M = 0.5 and in Fig. 7c for M = 0.9. The progressive 
shortening of the wavelength upstream of the source with increasing M and the lengthening 
downstream of it are manifestations of the Doppler frequency shift that occurs in an Eulerian frame 
moving with the noise source and not with the flow. 

From the point of view of achieving our objectives, computations are needed for jets issuing 
into stagnant or at most low velocity air. Thus, the acoustic boundary conditions will involve small 
values of Mach number which present no difficulty to us even though M can be arbitrarily large in 
the interior of the computational domain. 
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The semi— Lagrangian approach was then extended to the three-dimensional case. Results 
using upwind differencing for this method at a Mach number, M, equal to 0.75 are shown in Fig. 8a. 
Note that the sound field upstream of the monopole source is quite damped here due to numerical 
viscosity, which disagrees with theory. The theoretical result for a source fixed in a moving stream 
is that the magnitude of the sound is the same in the upstream and downstream directions, but varies 
in other directions. In Fig. 8b we plot the three-dimensional result using the new semi-Lagrangian 
method for a Courant number, C A , equal to 0.5 and a frequency that produces a wavelength that is 
half the size of the distance between opposite walls of the computational domain in the M = 0 case. 
The same frequency is used for M = 0.75, but the wavelength will change with location in a mean 
flow. Some improvement is seen compared to the result of Fig. 8a, although the new result does 
contain some damping upstream of the source. Further improvements are obtained by operating 
either at half the frequency (double the wavelength) as seen in Fig. 8c or at half the value of C A as 
seen in Fig. 8d. 

The major difficulty with using this approach is that the wavelength keeps getting shorter for 
a fixed frequency as M approaches 1.0. In this limit the number of points per wavelength for 
upstream propagating waves is small and the numerical inaccuracies in both the finite difference 
integration scheme and the interpolation procedures become large. The only clear way to improve 
the accuracy is to use a finer grid, which is what effectively happens in Fig. 8c. Decreasing the 
Courant number improves the accuracy of the interpolation procedures for the moving fluid case, but 
major improvements are achieved by using a finer grid. 

From a practical point of view inaccuracies in the preceding results for sound waves traveling 
upstream against a jet flow may not be detrimental to the overall jet noise results. This is because 
in a real jet of finite thickness the boundary between the jet and the ambient medium will refract 
sound waves traveling upstream back toward the jet axis. In particular, waves traveling upstream at 
sufficiently shallow angles to the jet axis will be trapped. The only way that this sound can reach the 
far field is if it is reflected by the nozzle or scattered by turbulent eddies. 


D. ADAMS-BASHFORTH TIME INTEGRATION 

The previous analysis has been based on standard finite difference techniques for constant 
time step increments. However, since we must couple the acoustics routines to the NEKTON flow 
code, which operates in general with variable sized time steps, it was decided to change the acoustic 
formulation to also treat variable sized time steps. We use the Adams-Bashforth method, the same 
approach used for time integration for the NEKTON flow code. The only difference is that 
NEKTON uses a third order method and we will employ a second order method. Although the second 
order method is less accurate in principle, we found it to be quite good for the range of parameters 
of interest to us and much more robust. For example, the third order method produced large errors 
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if the Courant number became too large but the second order method remained stable. Since the time 
step size, and thus the Courant number, was being controlled by NEKTON, we could not afford to 
risk getting into an unstable regime with the third order method. 

We will use the same philosophy as in the previous section of solving the wave equation in a 
frame moving with the local flow and then transforming results back to a fixed Eulerian. 

The Adams-Bashforth method considers first order differential equations of the form 

3p = f (21) 

dt 


with the solution at time step n+1 given in terms of values at step n and the earlier time n-1 as 


p(n + l) = p(n) + af(n) + bf(n-l) 


( 22 ) 


One can independently write a Taylor series expansion for p as 


2 2 

P(n + 1) = p(n) + At-^p(n) + ^ L -^-p(n) 
o t 2 dt 2 


(23) 


where At is the time interval between steps n and n+1. Substituting for the derivatives using Eq. (21), 
this becomes 


p(n + l) = p(n) + Atf(n) + |^(n) 


(24) 


To find an expression for f(n-l) in Eq. (22) another Taylor series expansion yields 


f(n-l) - f,„, - **00 ^(n) 

dt 2 st 2 


(25) 


where Atj is the time interval between steps n-1 and n. Then substituting for this in Eq. (22) 


p(n+l) = p(n) + af(n) + b f(n)-At x 


|f(n)-A tl g(n)j 


(26) 
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Equating the coefficients of ((9f)/(9t))(n) between Eq. (24) and Eq. (26) yields 

b = - (27) 

2At 1 


and equating the coefficients of f(n) yields 

a = At (1 + At/2Atj) ( 2g ) 

where At is the current time step and At x is the time step for the previous time. In the special case 
when At = At 2 


a = 3At/2 (29) 

b = -At/2 

so that information at the more recent time receives greater weighting. 9 

Our second order differential equation must be rewritten to utilize the form of Eq. (21) as 

d 2 P . df _ h (30) 

at 2 3t 

so that 


dP _ f (31) 

at 


which is identical to Eq. (21). The quantity h in Eq. (30) contains both the Laplacian of p, which is 
to be expressed in standard second-order finite difference form as in Eq. (10) for constant spatial 
increments, and the noise source term. The order of solution is important and is solved in the 
sequence 


f(n+l) 

= f(n) +ah(n) + bh(n-l) 

(32) 

p(n+l) 

= p(n) + af(n) + bf(n-l) 

(33) 
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Note that the same values of a and b apply to both equations. In terms of the Lagrangian solution this 
presents the same degree of difficulty as for the constant time step case since information is required 
at the n-1 time step and must be found by proper interpolation between the fixed and moving frames. 
The details of the integration process are contained in subroutine intgw. 


E. ROUNDAR Y CONDITION WITH FLOW 

The simple boundary condition given by Eq. (14) must be modified when there is flow 
through the computational boundaries. The boundary conditions are applied in the fixed reference 
frame. The case treated here with only a z component of velocity could be generalized to include the 
other components. When the Adams-Bashforth method is employed the coupling to the boundaries 
is through the function h in Eq. (30). This function contains values of the pressure at the desired 
point plus immediately adjacent points. For this reason the value of h on the boundary also depends 
on pressure outside of the boundary. The pressure at points located at one grid point outside of the 
boundary are determined by a generalization of Eq. (14) for the case with a z component of velocity 

d _V + ( C + U N ) iE. = 0 (34) 

at N ax N 


where U N is the component of velocity normal to the boundary. In this way values of p outside the 
boundary can be updated at new times in terms of previous values at the same point and previous 
values on the boundary. The details of this procedure are given in subroutine bndw. 


IV. FAR FIELD EXTRAPOLATION 


The extrapolation of near field acoustic data to the far field has been previously discussed. 1 
Briefly, the procedure is to gather pressure data on a large plane at some sideline distance from the 
jet, perform Fourier transforms in both time and the two space coordinates in the plane, and utilize 
a Green’s function solution of the Helmholtz equation (the frequency decomposed form of the wave 
equation) to obtain the sound spectrum at more distant locations. In the present case the Green’s 
function for the zero pressure boundary plane has been chosen. The far field solution is then given 
in terms of a distribution of dipoles on the plane whose strength is related to the transforms of the 
pressure data. This approach, which was developed and tested during our Phase I work, is presented 
in Appendix A. 
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Our studies show that accurate extrapolation of multipole source noise can be achieved if the 
distance between the source and the extrapolation plane is at least one— Quarter of the acoustic 
wavelength. 

Furthermore, if the jet axis coincides with the z axis and the far field observer is in the y-z 
plane passing through x = 0, then only the k x = 0 component of the data collected on the extrapolation 
plane at constant y is needed. Other wavenumber components produce far field sound outside of the 
y-z plane. This greatly reduces the computations since obtaining this component of k x results in a 
simple summation of all the pressure data as a function of x at each value of z instead of performing 
a complete Fourier transform. Test cases for a monopole point source show that reasonably accurate 
results can be achieved by performing the summation (or integration) over an absolute distance in x 
equal to ± one acoustic wavelength from the source. 

Methods for increasing the angular range in the y-z plane for accurate far field extrapolation 
of near field acoustic data are needed. A standard data processing technique is to extend the size of 
the spatial window by filling in the missing data points with zeros. 10 In general, the larger the size 
of the data extrapolation surface the broader the far field angular range that can be covered. 
Approximating the data by zeros should have some validity at distances far from the source. 

The second approach is to obtain the field strength at the desired wavenumber by 
interpolating values between wavenumbers computed from an FFT. In general the largest FFT 
wavenumber component, k p , parallel to the extrapolation surface that is still smaller than the total 
acoustic propagation vector, k, defines the largest angle of propagation from the normal to the 
surface. The wavenumber component, k N , normal to the extrapolation surface is given by 

k N = (k 2 - k 2 ) 1 / 2 (35) 

and must be real for propagating waves. The far field propagation angle, a, relative to the normal 
to the extrapolation plane is 


a = sin 1 (kp/k) (36) 

Larger FFT values of k p exist but they do not match conditions for far field propagation. 
Values of the Fourier transform for wavenumbers lying between k and the largest FFT value of k p 
that allows propagation are obtained by performing a linear interpolation between this wavenumber 
and the next largest FFT component that does not allow propagation. 

In Fig. 9a we show the ratio of the extrapolated acoustic intensity to the theoretical value for 
a monopole point source located at a normal distance from the extrapolation plane of 1/20 of its 
length and width and oscillating at a frequency that produces an acoustic wavelength that is 1/10 of 
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its length and width. The solid curve uses the Bingham window 10 in space and time. Rectangular 
window results are similar but display larger variations. The points to the left of the vertical line are 
obtained by the standard method and those to the right are based on interpolation between adjacent 
propagating and nonpropagating FFT components. Note that the first point to the right of the vertical 
line is quite large, indicating that more care is needed when examining wavenumbers near the 
propagation cutoff limit. The dashed line represents the results with zeros added to double the length 
and width of the extrapolation plane. Good results are obtained out to 64.2° from the normal for the 
standard case and out to 71.8° for the case of adding zeros without using the interpolation technique. 
Results for a reduced frequency producing a wavelength that is 1/5 of the length and width of the 
extrapolation plane display a narrower angular range as shown in Fig. 9b. The corresponding angles 
are 53.1° and 64.2°, respectively for the standard and added zeros cases. 

Our concern is in being able to accurately compute low frequency jet noise which is expected 
to have a significant contribution at small angles to the jet axis (large angles to the normal to the 
extrapolation plane). A range of 70° from the normal (20° from the jet axis) should be quite adequate 
for jet noise studies. To put this in perspective the wavelength. A, for jet noise can be obtained from 
the frequency expression 


and 


as 


fD/U = St 


(37) 


A = c/f 


(38) 


A 


D 

NfSt 


(39) 


where St is Strouhal number, U is the jet speed, c is the speed of sound, M is Mach number, and D 
is the jet nozzle diameter. From experiments the peak jet noise frequency occurs for St around 0.2. 
Arbitrarily picking M = 1 and the low frequency range to extend down to St = 0.1, then A = 10D. 
Thus an extrapolation plane would have to be 100D plus the distance from the jet exit to the axial 
location of the low frequency sources in length to get a ratio of 10:1 for the extrapolation dimension 
to wavelength ratio described in Fig. 9. Improving our extrapolation technique for low frequencies 
would enable us to obtain good results at useful angles with a smaller acoustic computational domain. 
An alternative would be to perform a separate integration of the wave equation over a coarser grid 
for low frequencies. 
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The extrapolation procedure described above is contained in a post-processing program called 
extrpfft which reads in data generated by the subroutine intgw in the main solution program and 
interpolates it onto an equally spaced time grid for each value of the axial coordinate z. An efficient 
way to run this code is to choose the number of interpolation time steps to be a power of two to 
facilitate Fourier transformation in time. Choosing the number of time steps in the acoustic 
integration procedure properly can also help. For example, in our runs the NEKTON integration time 
step did not change in size. On Voyager we were able to perform 320 NEKTON integration steps in 
each run, but not much more. The acoustic calculations were done every fifth time step so that we 
obtained results at 64 time points. Thus when running extrpfft we chose to obtain a 64 point FFT 
that yields 32 discrete frequencies. This spacing in time stepping is still small enough to satisfy C A 
< 0.2 to provide accurate results. 

In the present case only the lowest frequencies are of interest for jet noise. What we really 
need are sufficiently long data records to achieve a small enough frequency resolution to resolve the 
relatively low frequency components of jet noise. In the above case the 320 NEKTON steps 
correspond to 4.8 time units and the frequency resolution is then the inverse of this number. In terms 
of Strouhal number with velocity being one unit and diameter two units, the frequency resolution is 
St = 0.4167. This is not fine enough for jet noise so that a longer data record is needed. This can be 
readily achieved by stringing together the data outputs of a number of consecutive runs which are 
restarted with initial conditions equal to the final conditions of the preceding run so that the data is 
continuous. For example, combining four individual records into a single record and running the 
extrpfft program allows one to resolve Strouhal numbers as low as 0.1042. The combined time for 
four records corresponds to the time that it takes for the flow to travel ten jet diameters, or about two 
potential core lengths, at the inflow jet velocity. Stringing these files together requires that the very 
first line of data, containing run parameter information, be eliminated for all but the first data set. 
This combination can be accomplished because each additional run of the combined flow and 
acoustics codes is started with the flow and acoustics data of the last time step of the preceding run. 

Building up statistical reliability in the spectrum requires collecting a large number of 
independent data records of sufficient length in time. For this reason a lot of computer time is 
needed to obtain accurate jet noise calculations. Thus, to obtain ten records long enough to resolve 
St = 0.1042 requires 40 individual runs on Voyager. One way to reduce the calendar time is to submit 
a number of separate jet flow cases which do not have any of their time histories in common with any 
other case. This would be similar to performing laboratory experiments on a number of different jet 
nozzles having the same geometry and mean inflow conditions and averaging the data in the sense of 
an ensemble average. This can be accomplished by saving old flow restart files that cover the 
complete time needed if only one jet flow case were used or starting each flow case with some 
different random variation in the inflow velocity. 
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One must select the distance of the extrapolation plane from the center of the coordinate 
system in intgw prior to performing any runs. The computer guide describes how this is done. The 
lower the frequency of interest the further the extrapolation plane must be from the jet centerline. 


V. SELECTION OF FLOW PARAMETERS AND SPATIAL AND TIME SCALES 


Using the definitions of the acoustic Courant number and the Mach number one can write 


M 


” a *a 

c a Ax a 


(40) 


where both At A and Ax a are constants in the wave equation algorithm. By setting the maximum 
value of u, M is changed by varying c. We choose Ax a to be the largest scale that can resolve 
important small scale features of the turbulent noise source term to minimize the total number of 
computational points. Then high Mach numbers are achieved by decreasing C A or increasing At A . 
However, increasing At A too much will cut off our high frequency resolution. Therefore, a good 
choice is to lower C A . The distance, d, traveled by a sound wave is on the order 

d ~ cNAt A (41) 

d ~ C a NAx a 

where N is the number of time steps. Although d becomes smaller as C A decreases to raise M, the 
value of d needed to perform extrapolations of near field pressure to the far field also decreases. 
Thus, if the value of d needed to perform a good extrapolation scales with wavelength, 1 which is 
proportional to c, then a smaller value of d is not a problem. 

To achieve low values of M one would generally lower At A in Eq. (40) since C A has an upper 
limit. However, we do not recommend solving the differential equation form of the wave equation 
for very low M since we would question the numerical accuracy of the pressure calculations. This 
is because cancellations between different portions of quadrupoles become important and small 
numerical errors will dominate. In this case the integral solution approach would be preferable. 

Taking the initial jet radius as a unit length, the computational flow domain is 40 units long 
in the mean flow direction z and 20 units in width and height in the x, y directions. For z > 20 units 
from the inflow boundary the viscosity increases exponentially, reaching a value at the outflow 
boundary 1000 times larger than at the inflow. This forms a sponge layer that damps out disturbances 
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at the outflow boundary that could produce reverse flows that produce instabilities. The 
characteristic time is the unit length divided by the characteristic jet speed which is chosen to have 
a value of one. The velocity falls in the radial direction to a uniform value of 0.05. Instantaneous 
values of the source term s are computed using differentiating routines in NEKTON, and these values 
are interpolated onto the uniformly spaced grid used for the acoustic computations. The acoustic 
computations use the second order Adams-Bashforth time differencing scheme described previously 
to keep in step with variable NEKTON time steps. 

The acoustic computational domain has the same x, y dimensions as the flow domain, but the 
z dimension is somewhat longer as seen in Fig. 10. The grid spacing is 1/2 of a flow unit length (one 
jet radius) and the same in all directions. The z domain is divided into 128 points, corresponding to 
a length of 63.5 in z. Since only 81 points are needed to cover the NEKTON flow domain in z, the 
remainder are used as upstream and downstream overhangs to increase the length of the extrapolation 
plane. 


Thus, for M = 1 and St = 0.2, A/D = 5. To see where we stand with respect to the extapolation 
example of Fig. 9, we find with D equal to 2 length units, that A = 10, and 10 A = 100. Thus, the 
acoustic domain, which is 63.5 units long, falls somewhere between the extrapolation plane cases 
discussed in the last section for values of L/A equal to 5 and 10. 

The extrapolation plane is located at y = 5, corresponding to a distance of A/4 from the jet 
centerline when St = 0. 1 . The side walls of the domain are thus located at a distance of A/2 from the 
centerline at this Strouhal number. The ratio of the wavelength to these distances doubles for St = 
0.2 so that y can be reduced if desired. 

On the high frequency end, limitations are imposed by the resolution of the grid for the 
solution of the wave equation. Using the second order Adams-Bashforth technique a resolution that 
provides ten grid points per wavelength yields good results, and we have obtained reasonable results 
with six points per wavelength. However, dissipation becomes excessive for shorter wavelengths. 
This is not in itself a bad thing since it means that the computational technique will automatically 
rolloff the amplitude of high frequency, short wavelength, results which effectively provides an 
antialiasing filter for Fourier frequency analysis. For the present computations, with a grid spacing 
of two points per NEKTON length unit or four points per jet diameter, a high frequency criterion 
of six points per wavelength corresponds to A/D =1.5 and St = 0.667. This is in the range of the peak 
of the noise spectrum at 90° to the jet axis. Higher frequencies can be obtained with finer grid 
resolutions. Our resolution is presently not limited by the computation time required to solve the 
convected wave equation, but by the time needed to compute and interpolate the noise source term 
onto the acoustic grid. Expected increases in the code speed will lead to an ability to analyze higher 
frequencies. 
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A grid spacing of 1/2 provides adequate resolution of s for the larger values of z but not for 
smaller values where the derivatives of velocity that make up s and their spatial distribution can vary 
significantly. Since we are probably not interested in the smallest scales, a scheme was devised to 
smooth out these variations by performing a time average of s at each point. This is essentially the 
same as a space average since flow quantities are convected past a fixed point. We perform an average 
over 20 time steps, each with a time increment of about 0.015 time units so that the effective Strouhal 
number for updating s is 6.0, much higher than any frequency of interest for jet noise. The wave 
program is run every fifth NEKTON time step to keep the Courant number less than 0.2 

The extrapolation routine provides the pressure spectral level as a function of far field angle 
from the jet axis at each frequency. Polar (constant radius) results as well as constant sideline 
distance results are printed out. Since the angles do not repeat themselves (except for 90°) at different 
frequencies, a typical spectral plot at some designated angle must be obtained by interpolating plots 
of level vs. angle at a number of fixed frequencies and then replotting level vs. frequency at constant 
angle. The pressure is in a nondimensional form so that jets of different sizes have the same spectrum 
level and spectral content when frequency is expressed in terms of Strouhal number and distance from 
the nozzle in nozzle diameters. 


VI. RESULTS 


The flow studied was that of a round incompressible jet at a Reynolds number of 10,000 based 
on the diameter of the constant velocity portion of the upstream flow boundary. The spectral element 
mesh is in the form of a spider web design containing 13 elements in each plane perpendicular to the 
z axis as shown in Fig. 3. Element boundaries along the jet axis are atz = 0, 1,3, 5,9, 14, 21, 31, and 
40. 


The mean velocity at the inflow boundary, z = 0, was chosen as 

U, = 1.0; r = (x 2 + y 2 ) 1/2 < 1 


(42) 


and 


U, = 0.05 + 0.95 exp(-(r - 1.0) 2 ) (43) 

otherwise. The mean x and y components are taken as zero on the inflow and sidewall boundaries and 
U = 0.05 on the sidewalls. Thus, the inflow jet has a flat profile out to r=l and falls to an ambient 
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value of 0.05 at large distances. The Reynolds number is, with velocity and jet diameter set from 
above, determined by the specification of the viscosity. 

A random velocity in space and time was superimposed at the inflow plane on the mean 
velocity profile to help trigger flow instabilities that would lead to the development of turbulence. 
Flow computations showed that a very large turbulence level was produced. 

Next, a somewhat more coherent fluctuation of the radial components of velocity at the inflow 
plane in space and time was superimposed to try to break up the large scale jet structure. The spatial 
structure consisted of azimuthal modes of order 0, 1,2, and 3 and the time variation covers the range 
of frequencies 0 < St < 2. Both the mode and frequency structure changed slowly in a random manner 
based on the time scale of the basic fluctuation. This inflow condition is specified in subroutine 
ranin. This helped a little, but the turbulence level was still too high. 

Finally, fluctuating vorticity was added to the shear layer near the inflow boundary through 
the following conditions on the fluctuating z component, u B 

u B = 0.3 sin(Tr(r-1.707)/0.3) cos(2*t); |r - 1.707| < 0.3 (44) 

and u B = 0 otherwise. The above relation describes a shear flow that is in the form of a sine wave in 
r centered at the point of maximum mean shear at r = 1.707, having positive values on one side of this 
point and negative values on the other side. The phase, and thus the sign of the velocity perturbation, 
changes with time at a Strouhal number of 2. The disturbance looks like a ring vortex whose vorticity 
changes in magnitude and sign with time. 

The idea was that the vortical perturbations would be better suited to coupling with the flow. 
This appears to be the case since the general turbulence level has been reduced to realistic values on 
the order of 15% of the inflow jet velocity. Mean velocity profiles are shown in Fig. 11, turbulence 
profiles in Fig. 12, and the mean centerline velocity in Fig. 13. Note that the potential core is a little 
longer than might be expected. This might be due to the low Reynolds number or the fact that the 
initial jet shear layers are fairly thick. Typical curves of the value of the source term used to force 
the wave equation at one instant in time are shown in Fig. 14. 

The structure of s appears to be more like a number of traveling waves rather than a series 
of point sources. For that reason we speculate that our criterion for locating the extrapolation plane 
might be overly restrictive. In the case of point sources the wave number spectrum is very broad 
because of the small size of the source, and thus the range of phase velocities at constant frequency 
is also very large. Since most of this range corresponds to subsonic phase velocities whose energy 
predominates in the near field, greater accuracy for acoustic computations is obtained by collecting 
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data for extrapolation at larger distances. However, if the phase velocity of the turbulent noise source 
terms is already a significant fraction of the speed of sound, than a larger percentage of the energy 
will correspond to propagating modes. In that case data collection could occur closer to the source 
region without sacrificing accuracy. This is most true for far field propagation angles close to the jet 
axis where acoustic phase velocities projected parallel to the jet axis are close to sonic speed and least 
true for propagation at 90° to the jet axis where the required phase velocities approach infinity. 

Shown in Fig. 15 is a series of plots over a range of successive times of the k x = 0 component 
of the pressure on the extrapolation plane. This data is used directly to obtain the far field spectrum. 
At present we believe that these values are too large. This is probably due to values of the source 
term that are also excessive. Since the source term depends on differentiation of the instantaneous 
velocity field, more consistent values of the source term, for the frequency range of interest, can be 
obtained by filtering out small scale spatial variations prior to differentiation. This should be a top 
priority in any future work. 

A far field spectrum (Fig. 16) was obtained at 90° to the jet axis by combining data from four 
consecutive sets of near field pressure data into one long data record and applying the extrapolation 
routine. The spectrum peaks here at too low a frequency and at too high a level. This is a reflection 
of the spatial distribution of the source function in Fig. 14 which increases in amplitude with z. From 
simple dimensional arguments for a high Reynolds number jet we would expect that s should decrease 
with z. This is because the maximum turbulence level should change little with z, but the length scale 
increases with z so that velocity derivatives should decrease with increasing z. The behavior observed 
here is more typical of a low Reynolds number jet where the turbulence builds slowly with z and then 
bursts into a large amplitude fluctuation. Moreover, this growth in turbulence strength is seen in Fig. 
12. Methods to remedy this problem will be discussed shortly. 


VII. ADDITIONAL FEATURES 


While most of our attention has been directed at finding the pressure spectrum we note that 
other information is available. A postprocessing routine called read3 can be used to plot pressure 
field data generated by intgw. One can look at pressure along a specified line, e.g., constant x or z 
in the y = 0 plane, at a series of times or obtain a quasi three-dimensional plot of the pressure field 
in the y = 0 plane at one instant in time by importing the data into a spreadsheet program. 

The NEKTON post-processer enables one to plot a variety of flow information at instants in 
time on surfaces and along arbitrary lines as well as obtain complete time histories at selected points. 
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VIII. RECOMMENDATIONS 


A. COMPUTATION SPEED 

A number of improvements could greatly increase the accuracy and utilty of the computer 
programs. At the top of the list is reduction of computer time since this can be used to increase the 
spatial resolution, Reynolds number, and general computation throughput. We believe that this is 
possible for both the flow and acoustics codes and their interconnection. A higher speed program will 
be delivered after the official end of the contract period. 

Progress in speeding up the flow code is anticipated by modifying the method of solving the 
incompressible pressure equation in NEKTON which is in the form of a Poisson equation. It is 
presently solved using an iterative technique which consumes a considerable amount of time. Time 
reductions are possible by performing a direct computation, which requires a great deal of memory. 
We believe that Voyager should have enough available memory to facilitate this approach. 

We have found that replacing one routine in NEKTON with a library routine on Voyager can 
also speed up both the flow code and the interface between the flow and acoustics codes that 
computes the noise source strength. Thus, prospects are good that the estimates made in Section IV 
for the number of runs needed to produce a noise spectrum will be greatly reduced because many 
more time steps will be obtained in any single run. 

One weak point in the presented data analysis is the relatively course grid used for 
interpolation of the noise source function. An increase in computer speed would enable us to use a 
much finer grid, and thus provide a more accurate representation of this function, which is a sum of 
a product of derivatives. This is not so much of a problem if one is only interested in the larger 
turbulence scales responsible for low frequency noise. However, one must be careful that small scale 
fluctuations do not contaminate the low frequency results. This is essentially an aliasing problem in 
space and wavenumber rather than one in time and frequency. This can be alleviated by smoothing 
the flow velocities in space prior to computing the derivatives that are used in s. This smoothing is 
currently done for one velocity component in subroutine fmachw, and the same approach could be 
utilized for s. 


B. COMPUTER MEMORY 

Using a much finer grid in the acoustic computation can create other problems in terms of 
required memory since our finite difference solution uses an equally spaced mesh throughout the 
entire computation domain, e.g., halving the mesh spacing in a three-dimensional problem increases 
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memory requirements by a factor of eight. In reality a coarser mesh could be used outside of the jet 
region, but current thinking among researchers is that changes in spacing cause reflections of acoustic 
waves. We are not convinced that this is an ultimate barrier, but we think that there may be a more 
expedient approach that can be used for the jet noise problem. 

In our far field extrapolation technique only the k x = 0 component of pressure is needed on 
the extrapolation plane. A logical question would then be if this component can be solved for directly 
in the acoustic differential equation. The answer appears to be that it can if certain assumptions can 
be made. 

One proceeds by Fourier transforming the entire Phillips equation (or Lilley equation if 
desired) in x to obtain its spectral representation in k x . In the case when the x dependence of all the 
jet velocity components can be neglected then the only x dependence occurs in pressure. Fourier 
transformation in x is then a simple linear process with all x derivatives being represented by factors 
of k x . These x derivative terms then vanish if we are only interested in k x = 0. This leaves a partial 
differential equation in y, z, and t. Solution on a two-dimensional grid in y and z then reduces the 
memory requirments by orders of magnitude compared to the three-dimensional grid. 

In the general case velocity will be a function of x. We then recommend breaking up the x 
domain into a series of sections in which velocity is constant in x, but variable in y, z, t. One then 
numerically solves the partial differential equation in y, z, t for the k x = 0 component of pressure in 
each sector and sums the results of each sector. The size of each sector can probably be larger than 
the grid spacing used in the acoustic integration routine and these sectors need not extend out any 
larger in x than is needed to contain the noise sources. 

Physically, the assumption of no x variation in velocity in a sector is that the effect of 
scattering by the turbulence on that wavenumber component due to an x dependence is negligible. 
However, scattering by velocity variations in y and z is still included. In our current analysis we have 
emphasized formulations that focus on the effects of convection and refraction of sound by the flow 
which we believe to be dominated by the z velocity component. We think that this is a reasonable 
approach to try since scattering works in two ways: waves directed in one direction of interest can be 
scattered into other directions and waves travelling in other directions can be scattered towards the 
direction of interest. For a round jet with an axially symmetric average sound field these two 
scattering processes should balance. 


B. NOZZLE GEOMETRY 

Our original motivation for setting up the three-dimensional jet inflow conditions on a cross 
plane with no material nozzle boundaries present was that the jet turbulence a short distance 
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downstream of a nozzle would be relatively independent of its geometry for a high Reynolds number 
flow. In essence, any upstream time-dependent fluctuation would quickly trigger shear layer 
instability modes that would grow into turbulence and forget its origin. However, for the Re = 10,000 
jet studied the turbulence was strongly dependent on the nature of the upstream fluctuation. This 
dependence may in part be due to the low Reynolds number and the relatively thick initial shear layer 
width chosen to allow us to resolve the eddy structure. The growth of instabilities in thick shear 
layers will be susceptible to a smaller range of disturbances in frequency than for a thin layer. 

Two-dimensional slot nozzle studies which incorporated finite thickness nozzle lips showed 
the rapid development of a turbulence structure. For these reasons we strongly recommend that 
three-dimensional runs be performed with a nozzle lip. It is possible to insert a zero thickness lip 
which will develop a finite thickness boundary layer and will shed vorticity in response to velocity 
fluctuations occurring at downstream locations. 


IX. CONCLUSIONS 


The current status of numerical solutions of the Navier-Stokes equations is that the 
computation of time-dependent turbulent jet flows can be applied to current problems of interest. 
The major accomplishments are the development of codes for computing turbulent jet flows and 
acoustic propagation in time-dependent, nonuniform moving media and the connection of these two 
codes to form a jet noise prediction method. The present work constitutes a first step in the 
development of this approach. We believe that further refinements in the numerical code can lead 
to significant reductions in computational time as well as improved accuracy. The introduction of 
small, relatively inexpensive, parallel architecture computers will allow such computations to become 
more common. From this viewpoint anticipated hardware developments will revolutionize the utility 
of large scale numerical computations. 


28 



FIGURE 1 APPROACHES TO JET NOISE COMPUTATION 
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a. Equal Fin Lengths, Re = 300 



b. Equal Fin Lengths; Subgrid Turbulence Model 
Re «= 10,000 

FIGURE 2 VECTOR VELOCITY PLOT 







-10 

a. Grid in Cross-Sectional Plane 



b. Three-Dimensional Elements 
0 < z < 40 


FIGURE 3 SPECTRAL ELEMENTS FOR ROUND JET 
Units Referenced to One Jet Radius 
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DISTANCE FROM SOURCE IN GRID POINTS 

FIGURE 5 ACOUSTIC PRESSURE FIELD ALONG A LINE THROUGH MONOPOLE SOURCE 

Time Steps 1-40 
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-10 -8 -6 -4 -2 0 2 4 6 8 10 

DISTANCE FROM MIDPOINT OF SOURCE PAIR IN GRID POINTS 

FIGURE 6 ACOUSTIC PRESSURE FIELD ALONG A LINE THROUGH SOURCES OF DIPOLE 

Time Steps 100-200 
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WAVE AMPLITUDE 



x 


FIGURE 7 EFFECT OF MACH NUMBER, M, ON ONE-DIMENSIONAL SOUND WAVES 

Point Source Located at x = 0. 
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9 1-3A 




I I I I 

90 70 50 30 10 


ANGLE FROM JET AXIS, 0 
b. L/X = 5 

FIGURE 9 RATIO x OF PREDICTED TO THEORETICAL ACOUSTIC 
INTENSITY FOR POINT SOURCE 
L = Extrapolation Plane Length. 
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ACOUSTIC COMPUTATIONS 



FIGURE 10 RELATIVE DIMENSIONS OF COMPUTATIONAL REGIONS 
Units Referenced to One Jet Radius. 
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FIGURE 1 1 MEAN VELOCITY PROFILES 
R = Initial Jet Radius; D = Initial Jet Diameter. 
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U(r = °)/U jet RMS TURBULENCE LEVEL 
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TRANSVERSE DISTANCE, X/R 

FIGURE 12 TURBULENCE LEVEL 
y = o 
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FIGURE 13 MEAN CENTERLINE VELOCITY 
D = Initial Jet Diameter. 


38 




NOISE SOURCE STRENGTH 
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FIGURE 14 NOISE SOURCE STRENGTH 
y = 0; R = Initial Jet Radius; D = Initial Jet Diameter. 



FIGURE 15 PRESSURE ON EXTRAPOLATION PLANE 
k x = 0 Component; y/R = 5. 
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FIGURE 16 FAR FIELD NOISE SPECTRAL SHAPE 
90° to Jet Axis. 
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APPENDIX A 

NEAR TO FAR FIELD EXTRAPOLATIONS 


Obtaining far field acoustic results by increasing the size of the computational domain of the 
compressible flow equations would be exceedingly expensive. A more cost effective approach is to 
extrapolate near field results to the far field using asymptotic methods. The method of extrapolation 
we developed is based on the Green’s function solution of the Helmholtz equation, the time periodic 
form of the wave equation. The near field pressure can be computed either in terms of a frequency 
spectrum as a function of spatial location, a frequency-wavenumber spectrum, or space-time 
information. We will eventually require a frequency-wavenumber spectrum which, for example, is 
easily obtained from space-time data by the application of fast Fourier transform, FFT, techniques. 
These can be obtained in a matter of minutes on a personal computer, PC. This time is insignificant 
compared to the hours of computational time on a CRAY often required to obtain some compressible 
flow solutions. 

The Helmholtz equation for pressure, p, is 

v 2 p + k 2 p 
k 


0 

w/c 


(Al) 


where w is radian frequency and c is the speed of sound in a motionless, uniform medium. 

A relationship between the pressure on surfaces in contact with the medium and the pressure 
field is given by 11 


P 


J- f IgEE - p ^1 

4ir J ^ dn dn J 


dA 


(A2) 


where the integration occurs over a surface (or surfaces) denoted by A, dA is an element of surface 
area, d/dn signifies derivatives in the direction of the outward normal to the surface, and G is the so- 
called Green’s function satisfying the equation 11 


v 2 G + k 2 G = - 4x 6 (x - x 0 ) 


(A3) 
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The delta function on the right hand side is zero when the three-dimensional vector coordinate x # 
x 0 and has an integrable singularity in the limit x — ♦ x 0 . One solution to the equation is the free-field 
Green’s function G F given by 11 


G f 


exp(ik |x - Xpl) 


(A4) 


This is the solution for a point source of unit strength located at x 0 , and radiating into free space. 
However, there are an unlimited number of mathematically valid Green s functions obtained by 
adding any solution of the homogeneous Helmholtz equation. The importance of this is that there 
may be a Green’s function that will greatly simplify the evaluation of the integral equation for p in 
Eq. (A2). 

That integral contains values of both p and its normal derivative. The task would be 
simplified if G = 0 on the surface for two reasons. First, one of the integrations would be eliminated, 
but, more importantly, information on dp/dn from the numerical flow solution would not be 
necessary. This greatly reduces complexity and numerical errors. 

A case for which G = 0, on a particular surface is examined first and then generalized to 
arbitrary surfaces. Consider an infinite imaginary plane with one source located at position 

*0 = % + n 0®n (A5) 


in the half space above the plane where ? 0 is a position vector parallel to the plane, e n is a unit vector 
normal to the plane, and n 0 is the source’s distance above the plane. A second source is placed 
symmetrically below the surface at position 

*1 ■ % - n o®n 

and pressure is measured at arbitrary locations 

x = ? + n e n 


(A6) 


(A7) 


The following solution 
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^ exp(ik|x - x 0 |) exp(ik|x - Xj|) ( A g) 

|x - x 0 | |x - xj 

satisfies the delta functions at x„ in Eq. (A3) through its first term, which is identical to Eq. (A4). 
The second term, the solution for a negative source at x lt satisfies the Helmholtz equation, but does 
not lead to a delta function singularity as long as G is evaluated on or above the plane. 

To evaluate the integral in Eq. (A2) we must find the normal derivative of G relative to the 
source variable on the surface. Differentiation of the two terms in Eq. (A8) with respect to n 0 results 
in two quantities having the same magnitude. But the difference in the sign in front of n 0 combines 
with the difference in sign in front of the original terms to yield a result which is exactly twice the 
derivative of the free field Green’s function. 

This can be generalized to any surface with smoothly varying curvatures since the surface can 
locally be represented as a plane in the limit n 0 — * 0. The final result is 


P 




(A9) 


We apply this to the case of a given frequency dependent pressure distribution p(x 0 , y 0 , w) on 
a plane specified by z 0 = constant, and we wish to find the pressure at any point (x, y, z) above that 
surface. 


With 


we find 


|x - X(j| = ((x-x 0 ) 2 + (y - y 0 ) 2 + (z - Zo) 2 ) 1/Z 

d _ d 
drT 0 dT 0 


3Gj 

~dn~n 


dGt 

dn n 


dZn 


= ik 


exp(ik|x - x 0 |) 
lx - x 0 | 

L_1 f 

k|x - x 0 |J ^ 


z - z 0 

|X - XqI 


exp(ik|x - x 0 |) 

|x - x 0 I 


(A 10) 


(All) 
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Splitting the exponential term as follows 

exp (ik|x - x 0 |) = exp[i(k x x + k y y + k,z) - i(k x x 0 + k y y 0 + k^z,,)] 


(A 12) 


arbitrarily setting z 0 = 0 and defining the far field observation point as 


R = (x I 2 + y 2 + z 2 ) 1/2 
R ► x 0 , y 0 


(A 13) 


we obtain using Eqs. (A9), (All), (A 12) and (A 13) 


p(x,y,z,w) 


it fi . ±1 i 

2iR I kRj R 


I exp (ikR) 


(A14) 


L X L y 

1 = | | P (x 0 ,y D ,w) exp (-ik x x 0 - ik y y 0 ) dx 0 dy 0 


(A15) 


I = L x L y p k (k x , k y , w) 


(A16) 


The quantity I is the product of the effective area of integration and the two-dimensional finite 
Fourier transform in x 0 , y 0 . In the limit L x , L y -» oo, I is the infinite domain two-dimensional 
Fourier transform of the surface pressure. We recently found a similar approach by Veronesi and 
Maynard 12 who constructed Dirichlet and Neumann Green’s functions. Their Dirichlet Green’s 
function is the same as the expression in Eq. (A8). Note that our problem is not of the Dirichlet type 
since radiation boundary conditions are used at large distances from the source region. The acoustic 
intensity is proportional to the square of the pressure so that 


|p(x,y,z,w)| 2 


k 2 (L x L y ) 2 

4tt 2 R 2 



(|) 2 |p k (k x , k y ,w)| 2 


(A 17) 


A further simplification occurs when kR » 1, reducing the first bracketed quantity to unity. 
Using this and defining the pressure spectrum function P on the planar surface as 


p(k x ,k y ,w) = |I| 2 


(A18) 
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IP(X ’ y ’ Z ’" )|2 ■ ^5 ( I )' P ( k - k r “ ) (AI9) 

The appropriate far field values of k x and k y are obtained from geometric construction assuming 
sound travels along a ray from the origin x = y = z = 0to the far field point x F , y F , z F . Since 


k 


(K 



+ k 2 ) 1/2 = w/c 


(A20a) 


then 

K _ 

k ' (X* . yj * x*) 1 '* 

and 

ky = Vf 

k (x F + yf. + zp) 1/2 

In some applications it might not be practical to perform a single Fourier transform across the 
entire surface of interest, but it may be feasible to compute the transforms a number of times over 
smaller areas. If the pressure fluctuations are well correlated only over smaller elements of area, then 
the total area can be broken up into smaller portions and the contributions summed independently as 
follows: 


(A20b) 


(A20c) 


|p (x,y,z,w)| 2 


Pj (k x ,k y ,w) 


4jt 2 j 


R j 


(A21) 


Rt 


Note that Rj is the distance from the center of an individual element of area to the observation point. 
This approach enables one to compute the noise from surfaces that are arbitrarily large, but still 
satisfy the conditions Rj » l x , l y ; kRj » 1 where l x , l y are the extents of the area elements. Also, for 
a fixed number of points in an FFT algorithm one can treat higher wavenumber components of the 
pressure on the plane and a much larger range of far field directivity angles by analyzing a smaller 
region. 


We next tested the procedures just developed for extrapolating near field pressure data to the 
far field. The near field pressure was generated by combinations of simple point sources with known 
analytical solutions at any location. This near field pressure was Fourier transformed using FFT tech- 
niques, and the transform used to determine the pressure at some large distance using Eq. (A 19). 
This far field solution was then compared to the exact solution of the simple sources at the same 
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location. Calculations were performed on a PC compatible computer. A two-dimensional FFT 
algorithm obtained from a book by Newland 10 was utilized. All coding was done in FORTRAN 77 
and compiled on a MICROSOFT compiler. 

Consider first the imaginary plane defined by z = 0 and a simple point monopole source 
located at x = y = 0, z = -h having the form 


p = Re | exp (i (kR - Mt ))j 
k = cj/c 

R = (x 2 + y 2 + (z - h) 2 ) 1/2 


(A22) 


The complex quantity obtained by performing a Fourier time transform is 

p(«) = ex P < ikR) 


(A23) 


The two-dimensional transform in x and y is computed on the plane via our FFT routine with L x - 
L =1. When the far field observer is located at y F = 0, Eq. (A20c) shows that only the k y = 0 
component is needed. 

We define a nondimensional frequency parameter, W, which is the number of acoustic waves 
that fit along a segment of the ray path of length L x . Thus, the actual frequency in Hertz is given 
by 


f = 



(A24) 


For W = 100 we compared in Fig. A1 the ratio x of the predicted far field solution for the 
acoustic intensity to the exact solution for a simple source for z F = 100, y F = 0, h/L x = h = 0.6, over 
a range of x F values. The acoustic intensity is proportional to the square of the magnitude of 
pressure. Results for three different spatial data windows are shown. The rectangular window gives 
irregular results while the Hanning window^ gives excellent results for small values of x F , but trails 
off excessively for larger values. The Hanning window is zero at either end of the interval and rises 
to unity at the center. Its shape is that of a vertically displaced cosine curve that is always positive 
in value. 
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The Bingham window 10 gives the most desirable results since it covers almost as wide a range 
of angles as the rectangular window but is much more accurate than either of the two other windows 
over that range. It will be used exclusively in all the computations to follow in this report. Fourier 
transforms of 128 points each were calculated in the x and y directions to obtain these results. The 
Bingham window tapers the data in the first and last 10% of the interval with a cosine function. It 
is like a Hanning window with a central flat section in the middle that covers 80% of the interval. 

It makes sense when using the wavenumber representation on a plane that the ray connecting 
source to observer should pass through the plane, allowing a large enough area on it to describe the 
wavenumber structure. This is because the observer is viewing an image projected on the plane in 
our extrapolation technique, and the image will not be formed if the ray does not intersect the plane. 
In the latter test case the ray from the source which just hits the edge of the plane at y = 0 is at an 
angle of 39.8° from the vertical, where the computed pressure drops to zero. The Bingham window 
results differ from the exact solution by only 1.3 dB at 33°. To increase the useful angle we reduced 
h/L x to 0.1. We also had to decrease W to 50 since the 128 sampling points were not enough to satisfy 
the Nyquist criterion 13 when the rays were nearly parallel to the plane. In this latter situation there 
are not enough points to characterize the projection of waves on the surface when W = 100. In Fig. 
A2 we have a much improved angular spread with discrepancies of less than 1.5 dB out to 55°. 
Reducing h/L x to 0.05 yields a discrepancy of 1.36 dB at 55°, with large improvements relative to the 
h/L x = 0.1 case at larger angles. 

We also examined the computation of multipole sound fields. This is of interest because real 
sound sources often display such characteristics. The technical problem is that multipole sources have 
distinctly different near field and far field characteristics and the acoustic pressure waves of interest 
to us may be small in amplitude compared to nonpropagating hydrodynamic pressure fluctuations if 
one is too close to the multipole. For the monopole case we saw that it was beneficial to move a 
measurement surface of a given size as close to the source as possible to obtain good far field results 
over the widest angular domain. But there is a limit to how close the surface can be for a multipole 
since any numerical errors in the nonpropagating pressure will contaminate the acoustic portion of 
the pressure. 

To study this effect we compared the acoustic intensity computed from the Green’s function 
approach for dipoles against the exact solution. The dipoles consisted of two point monopoles 
oscillating sinusoidally, but 180° out of phase with each other. In one case, the perpendicular dipole, 
the line joining the two monopoles is perpendicular to the measurement plane, and in the second case, 
the parallel dipole, the line joining them is parallel to the plane and aligned with the x axis. The 
distance between the monopoles is kept small compared to the distance between any monopole and 
the measuring plane. 
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In Fig. A3 we see that for the perpendicular dipole case decreasing the value of h, the distance 
from the measurement plane to the center of a line connecting the two monopoles, generally improves 
the results over a wide range of angles up to about h/A = 0.24 where A is the acoustic wavelength. For 
h/A = 0.16 there is a small, less than 0.4 dB, but noticeable drop in the computed acoustic intensity. 
For a parallel dipole, departures from the exact solution are typically on the order of 1 db for h/A = 
0.24. Some of this error may be due to the fact that the exact solution includes both constructive and 
destructive interference effects between the two sources, causing rapid changes in the pressure’s 
spatial amplitude. Based on the results of Fig. A3 we recommend that the distance between the 
measurement plane and the source region be no less than about A/4. 

Based on these numerical experiments we can formulate a recommendation for the size and 
location of the measurement surface. Consider the source region to be distributed in the flow 
direction along a length L s . The near field measurement surface length should be at least 20 times 
the distance from source to surface for a single simple source, which means a length 10 times the 
distance upstream and the same length downstream of the source region for distributed sources. 
Defining the total length L T as 


L t = L + L s (A25a) 

then 

L t > 20h + L s (A25b) 

or, since h > A/4 

L T > 5A + L s (A25c) 

For the transverse case we would by analogy replace Lg by the width of the distributed source region 
to obtain the surface width. 

For an axially symmetric jet the source region will extend from r = 0 out to around r = D for 
a spreading flow. Thus, the measuring surface must be no closer than D + A/4 from the centerline. 
To meet the condition h > A/4 for points on the jet edge requires 

L x /D > 20 + 5A/D + L s /D (A26) 
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0 5 10 15 20 25 30 35 40 

ANGLE FROM NORMAL TO MEASUREMENT PLANE 

FIGURE A 1 EFFECT OF SPATIAL WINDOW ON FARFIELD 
PREDICTION OF MONOPOLE PRESSURE 

X, Ratio of Computed to Exact Acoustic Intensity; h. Distance of 
Source to Plane; L x , Length of Measurement Plane; h/L x = 0.6. 
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ANGLE FROM NORMAL TO MEASUREMENT PLANE 


FIGURE A 2 EFFECT OF MONOPOLE DISTANCE FROM MEASUREMENT 
PLANE ON FARFIELD PREDICTION 

X, Ratio of Computed to Exact Acoustic Intensity; 
h. Distance of Source to Plane; A, Wavelength; 

L x , Length of Measurement Plane; L x /A = 50. 
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ANGLE FROM NORMAL TO MEASUREMENT PLANE 

FIGURE A3 EFFECT OF DIPOLE DISTANCE FROM MEASUREMENT 
PLANE ON FARFIELD PREDICTION 

X, Ratio of Computed to Exact Acoustic Intensity; 
h. Distance of Source to Plane; A, Wavelength; 

L x , Length of Measurement Plane; L x /A = 16. 
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APPENDIX B 

COMPRESSIBLE TURBULENCE JET 


A. NUMERICAL SIMULATION OF COMPRESSIBLE TURBULENCE 
1. Dynamical Equations 

We have developed a compressible flow computer code that allows efficient simulation 
of a variety of low-to-moderate Mach number compressible flows, including jets and other free- 
shear flows. We consider the motion of a viscous, compressible fluid with a perfect-gas equation of 
state. The equations of conservations of mass, momentum, and total energy are 

d ± + _!_(p Ui ) = 0 (Bl) 

dt dx, 1 


3_ 

at 


(PUj) 


, _ <**.,> 


. _L ± fs u 

dx i Re 0 5xj [ J 



^ E T . (E T . p)u, - 


! v 2 t 

M*PrRe 0 (7 - D 


+ 


_2 a 

Re 0 3x 


A s « 



(B2) 


(B3) 


where p is the density, p the pressure, T the temperature, u = (u x , u 2 , u 3 ) the velocity, 

A = V»u 


(B4) 


the dilatation, and 


a Uj 

a U j 

+ 

axi 

dx j 


(B5) 
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the strain-rate tensor. All the field variables are non-dimensionalized in terms of a reference density, 
p 0 , velocity u 0 , length 1 0 and temperature T 0 . The bulk viscosity is assumed to be zero. We use the 
convention that repeated indices are summed over. Finally, the perfect-gas law is 


pT 
7M J 


(B6) 


The total energy density 


E x - + E k 


is the sum of the internal energy density 


El - 



and the kinetic energy density 


E 


K 


xM 2 


(B7) 


(B8) 


(B9) 


There are three non-dimensional parameters: the Reynolds number. 


- ^o u o l o (BIO) 

Ke 0 = 

the Mach number, which is defined by the ratio of a typical velocity u 0 and sound speed 
c 0 = hRTo) 1 / 2 


M 0 



(B 11) 


and the Prandtl number 
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P 


r 


C p" 


(B12) 


where 7 is the specific heat ratio, R the perfect-gas constant, Cp the constant pressure specific heat, 
fi the dynamic viscosity and k the thermal conductivity. 

The resulting non-dimensional equations of motion are 


dp 

at 


+ *T (pu i> 


0 




JLipu.^) - - 

♦ A(E t + p )u j 

a Xj 




dp + _2_ _d_ 
axj Re 0 ax k 

1 

M^PrRe 0 (7 

2 d (c 

Re 0 ax k I 



with the perfect-gas law 


pT 

'yMjj 


(B13) 


(B14) 


(B15) 


(B16) 


2. Initial and Boundary Conditions 

We solve Eqs. (B1)-(B3) numerically with inflow-outflow boundary conditions in a 
three-dimensional box geometry. A spectral method is used following the work of Spalart in which 
a polynomial subtraction is used after a coordinate mapping to represent the non-periodic flow as the 
sum of a low-order polynomial-like contribution and a periodic correction term. The periodic 
components of the dependent variables are then expanded in a spectral/Fourier series and spectral 
collocation methods are applied to solve the Navier-Stokes equations. 

There are several novel technical details in the above procedure. First, viscous sponge 
layers combined with parabolized equations of motion are used to implement efficient outflow 
boundary conditions. Variable viscosity in the outflow direction allows the computation of complex 
jet flows on finite-sized domains. 
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Second, a linearly varying tanh mapping is used to collapse the infinite jet flow 
geometry into a finite domain. This transformation gives an elegantly simple set of transformed 
Navier-Stokes equations that are efficiently implementable in the present spectral representation. 

Third, a sub-grid scale eddy viscosity transport model is included to enable large-eddy 
simulations of these flows. The sub-grid scale viscosity is that derived from renormalization group 
(RNG) methods. This RNG viscosity gives an effective treatment of the laminar-turbulent interface 
at the edge of the jet. In effect, the RNG model accounts properly for ’low’ Reynolds number 
effects, which occur where laminar fluid is entrained into the turbulent jet. 

Fourth, the initial conditions are chosen to represent a jet profile. This is typically 
done by either: (a) starting with a field of homogeneous compressible turbulence and adding in the 
mean jet flow field and ’shaping’ the turbulent fluctuations to have a jet-like intensity pattern; or (b) 
taking an incompressible jet flow run and using it as initial and inflow boundary conditions for the 
compressible jet run. As a more complete library of compressible jet flow simulations becomes 
available, these earlier runs can be used as the initial and inflow boundary conditions for compressible 
jet runs made under different operating conditions. 


B. DYNAMICS OF COMPRESSIBLE TURBULENCE 

The effects of compressibility become more and more significant in the dynamics of 
compressible turbulent flows as the fluctuating (rms) Mach number M increases. When the 
fluctuating turbulent velocity becomes comparable with the speed of sound, interactions between 
vorticity and sound waves are substantial. The most relevant parameter characterizing the 
incompressible turbulence is the micro-scale Reynolds number R v Turbulence of different kinds is 
conveniently compared with reference to R A . The statistical properties of turbulence, such as the 
form of the energy spectrum, the skewness and flatness of velocity-derivatives, for a given R A are 
more or less universal, irrespective of the kinds of turbulence. 

For compressible turbulence, there are additional parameters which characterize the property 
of the field. In our paper, 14 we studied numerically a forced compressible turbulence and showed 
that the statistical properties of the velocity field vary depending not only on the rms Mach number 
but also on the intensity of fluctuations of thermodynamic variables (e.g., density fluctuation Sp' / <p>) 
and the small- and large-scale compressive ratios r cs and r CL . Since the local Mach number may take 
much larger values than (typically twice as large as) the rms Mach number, the formation of shock 
waves is expected even at moderate rms Mach numbers (» 0.5). In this case, interactions between 
vorticity field and shock waves are important . 
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In Ref. 14, we decomposed the velocity field into the compressive and rotational components 
and examined the energy-transfers among these two components and the internal energy through 
various nonlinear interactions. The coupling between the compressive and rotational components is 
relatively weak compared with interactions among modes of the respective components. The spectrum 
of the rotational component of energy does not depend on M and obeys the same Kolmogorov 
similarity law as for the incompressible case. The compressive component, on the other hand, 
depends strongly on M and the spectrum has a longer tail for larger M. We observed a periodic 
energy-exchange between the internal energy and the compressive component of kinetic energy for 
high— compressive ratios. This energy-exchange occurs through the pressure-dilatation interaction 
term. A similar observation had been made by Passot and Pouquet 16 by their numerical simulation 
of two-dimensional decaying turbulence. The contribution to the energy-dissipation by the 
compressive component of velocity (longitudinal mode) is effective at high— Mach number and 
high-compressive ratios (see also Refs. 16 and 17). 

We also studied a role of shock waves in generation of vorticity. We found that vorticity 
production due to baroclinic generation is effective by curved shocks and that the pressure and 
density gradients are nearly parallel. This observation is different from the two-dimensional 
simulation by Passot and Pouquet. They observed that vorticity is produced by shock wave collisions 
near which pressure and density gradients are nearly orthogonal. The effects of compressiblity at low 
Mach number were analysed by Sarkar et al. 16 and Erlebacher. 18 They found that two different 
asymptotic states exist which depend on Mach number and pressure fluctuation of the initial velocity 
field. 


Here we study the statistical properties of a freely decaying compressible turbulence. We 
investigate the decay law of kinetic energy, the Kolmogorov similarity law of the kinetic energy 
spectrum, the relaxation of the compressive ratios to some universal values, and the energy-exchange 
among the internal energy and the compressive and rotational components of kinetic energy and so 
on. 


C. SUMMARY OF RESULTS AND DISCUSSION 

We studied some of the statistical properties of decaying compressible turbulence by solving 
numerically the compressible Navier-Stokes equation supplemented by the equation of state for an 
ideal gas. The initial flow field is a fully developed turbulent state which was generated by a random 
external force. The flow is confined in a periodic cube. The characteristics of the flow are different 
depending on the micro-scale Reynolds number R A , the rms Mach number M, the small- and 
large-scale compressive ratios r cs and r CL . We observed the following interesting results: 
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(1) The kinetic energy (and also the velocity squared) decays exponentially in time. This 
decay law is clearly seen for the case that the compressive ratios are so small that the energy-exchange 
between the kinetic and internal components through the pressure-dilatation interaction may be 
negligible compared with the energy-dissipation by viscosity. If the compressive ratios are large, the 
kinetic energy shows a decaying oscillation resulting from a periodic energy-exchange between the 
compressive component of kinetic energy and the internal energy. 16 The exponential decay law of 
energy has been observed in the incompressible turbulence too. 1 Thus, it is inferred that the 
exponential decay law of energy may be one of the universal properties of turbulence realized as at 
least a certain period of the decay process. 

(2) The kinetic energy spectrum (as well as the power spectrum of velocity) of the rotational 
component obeys the Kolomogorov similarity law unless the Mach number is not too large (M < 0.6). 
On the other hand, the spectrum of the compressive component varies very much at larger 
wavenumbers depending on the Mach number and the compressive ratios. This suggests that even 
if the Mach number is moderate (M « 0.5), the compressive and rotational components behave rather 
independently of each other as long as the modal distribution of energy is concerned. 

(3) In the freely decaying process of turbulence, the micro-scale Reynolds number and the 
rms Mach number decrease in time but the two compressive ratios seem to approach some universal 
values. Because of the shortness of the simulation period, however, it is not clear from the present 
numerical simulation whether the asymptotic state, if it exists, is unique or not (see Refs. 16 and 18). 

We have reported here some statistical properties of compressible turbulence with relatively 
small Reynolds number (R A < 40) and rms Mach number (M < 0.8). These values of parameters seem 
to be the upper limits reachable with reasonable accuracy in a simulation of 64 s collocation points. 
We have shown that the Kolmogorov similarity law does still hold for the rotational component of the 
flow. However, a bigger simulation is necessary to realize the Kolmogorov inertial range spectrum 
and to examine interactions among Fourier components for higher Reynolds and Mach number 
turbulence. 
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APPENDIX C 

COMPUTER OPERATION GUIDE 


In the description which follows, the jet flow code is referred to collectively as NEKTON 
whereas the acoustic code is referred to as WAVE. A brief synopsis of the important files appears 
below. See restrictions on the use of the code outside of NASA in Section G. 


A. BRIEF SYNOPSIS OF IMPORTANT FILES 

1. NEKTON Pre-processor 

The pre-processor is invoked by entering PRENEK. This is a graphical application requiring 
the use of a 2 or 3-button mouse. Available versions of PRENEK support SUN (Sparc), 
TEKTRONIX (4107), DEC(3100) and SILICON GRAPHICS (Irix) terminals. The pre-processor 
allows one to define the NEKTON grid and set initial and boundary conditions for the problem. The 
.inf, .cmn and .rea files mentioned below are also generated. For more complete information on the 
use and operation of the pre-processor refer to accompanying NEKTON manuals. 

a. NEKTON Pre-processor Input Files 

No input files are required by the pre-processor. 


b. NEKTON Pre-processor Output Files 

*.inf - contains information pertaining to inflow and viscosity profiles. File is renamed as 
a .f file as outlined below. 

*.cmn - contains variable definitions and declarations. Must be edited as outlined below. 

*.rea - contains time stepping information, NEKTON grid definition and various physical 
flow conditions. 
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2. NEKTON Post-processor 

The post-processor is invoked by entering POSTNEK. This is a graphical application 
requiring the use of a 2 or 3-button mouse. Available versions of POSTNEK support SUN (Sparc), 
TEKTRONIX (4107), DEC(3100) and SILICON GRAPHICS (Irix) terminals. The post-processor 
allows one to view spatial and time-history graphical representations of various NEKTON fluid flow 
parameters. The .rea, .fid and .his files mentioned below are required as inputs to the pre-processor. 
For more complete information on the use and operation of the post— processor refer to accompanying 
NEKTON manuals. 

a. NEKTON Post-processor Input Files 

*.rea - contains time stepping information, NEKTON grid definition and various physical 
flow conditions. 

*.fld - contains NEKTON field data of instantaneous velocity and pressure. 

*.his - contains NEKTON field data of velocity and pressure histories. 

b. NEKTON Post-orocessor Output Files 

No output files are generated by the post-processor. 


3. WAVE Post-processor 

Two (2) separate post-processing tasks are available. The one chosen depends on the 
information required. These post-processors are interactive and non-graphical. They require answers 
to a few self-explanatory prompts. 

a. WAVE Post-processor Input Files 

specs - contains tabulated information for pressure, source terms, etc. for entire X-Z plane. 
This file is used as input for the READ3 post-processor. 

ffts — contains pressure data at Y = YEXTRP as function of Z and time. This file is used as 
input for the EXTRPFFT post-processor to find the far field pressure spectrum. 
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b. WAVE Post-processor Output Files 

The output files are specified by the user in response to the prompts when the post-processors 
are invoked. 

The output from the READ3 post-processor is the data from "specs" along a constant-X or 
a constant-Z. Data is available to form a three-dimensional plot of pressure in the X-Z plane using 
standard graphical spreadsheet programs. 

The output from the EXTRPFFT post-processor lists tabulated values of the pressure 
spectrum vs. angle from the jet axis, at various fixed frequencies (Strouhal numbers) at both constant 
radius and constant sideline distance. 


4. NEKTON Flow Computation Files 

a. NEKTON Source Files 

try.rea - contains number of time steps, step interval, restart file name, element Cartesian 
coordinates, boundary types, flow simulation flags, number and location of history points, physical 
flow conditions such as density and viscosity. 

try.f - contains routines which describe more complicated inflow and viscosity profiles than 
allowed in preprocessor. 

try.cmn - contains definition and declaration of several NEKTON parameters and variables. 

ranin.f - contains definition of random boundary inputs. This is case specific and may be 
modified or excluded. Called from try.f. 

rancray.f - contains CRAY-2 random number generator calls. Called from ranin.f. 

mgdummy.f - contains dummy subroutine stubs for various NEKTON calls. 

nekton l.f - contains calls to acoustic routines WAVE and time-stepping routine calls. Main 
program DRIVE is located in this module. 

nekton2.f - contains interpolation and source smoothing schemes. 
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nekton3.f - contains form of fluid flow equation utilized by solvers (Stokes or Navier-Stokes) 
and consequent grid modelling. Contains viscosity corrections for certain flow conditions. 

Depends on the flags set in try.rea. 


b. NEKTON Input Files 

bounds.dat - contains values for bounds of NEKTON and interpolation domains. This file 
will be created if it does not exist. 

ielist.dat - contains values representing the NEKTON domain locations of the WAVE 
interpolation points. These values are calculated only once and saved for use in future runs unless 
new domains are defined. This file will be created if it does not exist. 

c. NEKTON Output Files 

source.dat - contains values for source terms at each point on the NEKTON grid. 

intsrc.dat - contains values for source terms at each point on the WAVE interpolation grid. 

intvel.dat - contains values for Z-velocity at each point on the WAVE interpolation grid. 

*.fld - contains NEKTON field data of instantaneous velocity and pressure. Used as input 
during post-processing. 

*.his - contains NEKTON field data of velocity and pressure histories. Used as input during 
post-processing. 

*.rst - contains more NEKTON field data consisting of mean and root-mean-square velocities. 

The *.fld and *.rst files must be renamed and submitted as restart NEKTON files for the 
immediately succeeding run. 


5. Acoustic Computation Files 
a. WAVE Source Files 

wavnek.f - formed by concatenation of FORTRAN source files wavw, intgw, initw, prntw, 
bndw, srcw, shftw, fmachw. Contains the acoustic solutions. Routines wavw, intgw and prntw called 
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from nekton2.f. Inputs are from NEKTON via common blocks specified in include file WAVC. 
Outputs files ffts, specs, prnt.rst and courant.chk. On restart prnt.rst must be renamed as init.rst. 
Init.rst now becomes an additional input for wavnek.f. 

b. WAVE Input Files 

init.rst - contains values for various parameters needed for WAVE restart. This file will not 
be present until after the first run since the user must rename the output file prnt.rst as init.rst. 

c. WAVE Output Files 

prnt.rst - contains values for various WAVE parameters at the time the current run was 
terminated. This file must be renamed as init.rst so that it will be used as input on the immediately 
succeeding run. 

specs - contains tabulated information for pressure, source terms etc. for entire X-Z plane. 
This file is used as input for the READ3 post-processor. 

ffts - contains pressure data at Y = YEXTRP as function of Z and time. This file is used as 
input for the EXTRPFFT post-processor to find the far field pressure spectrum. 

courant.chk - contains information pertaining to any time step for which the Courant number 
is 0.2 or greater. 


6. INCLUDE Files 

ALL - contains parameters setting the number of interpolation points, the first time step 
at which interpolation occurs and the number of steps over which time-averaged smoothing is 
performed. 

WAVC - contains parameters setting indexes for the acoustic domain, the interpolation 
region within this domain, and the Cartesian coordinates (referred to NEKTON space) of these 
boundaries. Also contains values for gamma and Courant number used in acoustic calculations. 

READ3C - contains parameters setting indexes for the acoustic domain. Indexes are the 
same as in INCLUDE file WAVC above. 

SPONGE - contains formula defining ’sponge layer’ and the viscosity distribution in this 

layer. 
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SIZE - contains parameters setting the dimensions of the NEKTON element grid and the 
order of the grid mesh resolution. 

TOTAL - contains parameters needed in every NEKTON routine. 

B. TO CHANGE RIJN PARAMETERS FOR GIVEN CASE 

Edit try.rea file as follows : 

NSTEPS - to change number of time steps 

IOSTEP - to change time step interval between acoustic calculations 

KINVIS - to change viscosity restart file - search for the string ’Reading’ and modify the 
filename on the same line on which the string occurs. 

IFMODEL - Should be set to T for all runs invoking the turbulent model. Should be set to 
F otherwise. 

IFKEPS - Should be set to T to obtain non-mixing length RNG solutions. Should be set to 
F otherwise. Value is ignored when IFMODEL is set to F. 

Edit try.f file as follows : 

SUBROUTINE INFLOW - to modify inflow boundary conditions. This requires knowledge 
of which elements and sides represent inflow boundaries. 

Edit try.cmn file as follows : 

LX, LY, LZ, LZ1 - to increase mesh resolution, that is, the number of nodes per NEKTON 
element. 

Edit SPONGE file as follows : 

Z,TAU1 - modify the z-value (in NEKTON Cartesian coordinates) at which ’sponge’ layer 
starts. Also change viscosity distribution in layer by modifying or rewriting formula. 
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Edit WAVC file as follows : 

MAXX, MAXY, XINT, XBND1 etc. - Modify dimensions of acoustic domain and NEKTON 
interpolation region. Parameter names are self-documenting: 

MAXX, MAXY, MAXZ - Overall dimensions of the WAVE domain 

XINT, YINT, ZINT - Overall dimensions of the NEKTON domain 

Note that the NEKTON domain dimensions are less than or equal to the WAVE 

domain dimensions. Note also that the NEKTON domain is centered within the 

WAVE domain whenever the two sets of dimensions are different. 

XBND1, XBND2, YBND1, YBND2, ZBND1, ZBND2 - Cartesian coordinates in NEKTON 
space of the bounds of the NEKTON domain 

XSTART, XEND, YSTART, YEND, ZSTART, ZEND - Cartesian coordinates in NEKTON 
space of the bounds of the NEKTON interpolation region. 

ACOUST - Change value to FALSE in order to suspend acoustic calculations. 

DUMPALL - Change value to TRUE in order to obtain NEKTON field dumps at every 
IOSTEP interval. These fields may be later post-processed if desired. A value of FALSE ensures 
only the final dump will be obtained. 

GAMMA, VELC - Change values for gamma (the ratio of specific heats) and Courant number 
used in acoustic calculations. The Courant number Cr is related to VELC by the formula Cr = VELC 
* DT/DZ where DT is the delta-time value passed to WAVE from NEKTON, DZ is the WAVE grid 
spacing, and VELC is the sound velocity referred to unit jet velocity. A change in any or all of 
VELC, DT or DZ must be such as to maintain Cr < 0.2. A change in DT is effected when the 
IOSTEP variable mentioned previously is changed. Changes in DZ occur when the domain dimensions 
and/or number of WAVE grid points are changed. VELC can be changed directly. 

FFTSTRT - Change the time at which far field extrapolation data starts being calculated. 

YEXTRP - Change the physical Y-coordinate of extrapolation plane. An increase in this 
parameter may improve the accuracy of the far field data. 

DIAJET, VELJET - Change jet diameter and velocity, respectively. These parameters are 
changed only if a new NEKTON grid and boundary specifications are defined. 
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ST - Change Strouhal number. This is really a maximum value estimate which can only be 
fixed from observation of the delta-time values. 

Edit EXTRAPIN as follows : 

MAXX, MAXY, MAXZ - Change parameters in line with changes made in WAVC. 

Edit TOTAL file as follows : 

VKINIT, VKTIME - Change parameters to adjust the viscosity at low Reynolds number, and 
the time at which the viscosity switches value for high Reynolds number flow. This technique 
produces large— eddy flow prior to the initiation of the RNG turbulence model to establish large scale 
velocity gradients. 


C. TO DEVELOP PARAMETERS FOR A NEW CASE 

1) Run NEKTON pre-processor by entering PRENEK. The instructions for running the 
pre— processor are found in the manual provided. 

2) Rename the .inf file as .f file. 

3) Edit the .cmn file as follows : 

Add the following line after line 3 of file : 

PARAMETER (MZ=LX1) 

4) Edit the .rea file as follows : 

Change value of IFMODEL to F 
Change value of IFKEPS to F 

5) Follow the procedures outlined previously under "To Change Run Parameters for Given 

Case". 

6) Submit job by using the makefile supplied. See details under "How to Submit a Job under 

UNIX". 
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7) At the end of each run save the .fid, .his, .rst files under new names, e.g., save 
try.fld,try.his,try.rst under try 1. fid, try 1. his , try 1. rst. The saved files may be used for NEKTON 
post-processing if desired. 

8) Inspect the first value on the first line of the .rst file. This value is the restart time. 

9) Edit .rea file as follows : 

If the restart time from instruction no. 8 above is equal to or greater than VKTIME set in 
include file TOTAL from instruction no. 5, then change values of IFMODEL and IFKEPS to T. 
Search for the string ’Reading’ and change the restart file name to reflect the names used from 
instruction no. 7. 

10) If ACOUST variable is set to TRUE, then acoustic results are available. Save files named 
ffts and specs under new names for WAVE post-processing if desired. Rename file prnt.rst as init.rst. 


11) Repeat instructions starting at instruction no. 6 for additional runs without modifications 
to the flow. Go to no. 5 to change basic parameters. Go to no. 1 to change geometry and repeat all 
steps. 


D. TO SET UP RUNS FOR EXISTING CASE 


Repeat instructions starting at instruction no. 6 as outlined previously in "To Develop 
Parameters for a New Case". 


E. TO POST-PROCESS DATA FROM NEKTON AND WAVE 

NEKTON Post-processing 

Three (3) files are needed by the NEKTON post-processor. These are the .fid, .his and .rea 
files. The post-processor is invoked by entering POSTNEK. For more complete information on the 
use and operation of the pre-processor refer to accompanying NEKTON manuals. 

WAVE Post-processing 

Two (2) separate post-processing tasks are available. The one chosen depends on the 
information required. These post-processors are interactive and non-graphical. They require answers 
to a few self-explanatory prompts. 
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Invoke one WAVE post-processor by entering READ3. The input file for the post-processor 
is the file named "specs". The output from the READ3 post-processor is the data from "specs" along 
a constant-X or a constant-Z for one value of Y. Data is available to form a three-dimensional plot 
of pressure in the X-Z plane using standard graphical spreadsheet programs. 

The Y value may be adjusted by changing the value of NCENT in INCLUDE file WAVC. 
The information is printed only for the last time step of each run. 

Invoke the other post-processor by entering EXTRPFFT. The input file for this 
post-processor comprises the modified and concatenated files ffts produced by each 
NEKTON/WAVE run as described previously. 

The "ffts" files are modified as follows : 

(1) Change the fifth field of the first line of the first file in the catenation chain such that 
it reflects the actual number of time steps represented when all the files are concatenated. It is 
strongly suggested that the number of files in the catenation chain be such that the actual number of 
time steps is a power of 2, i.e. 64, 128, 256, 512, etc. This results in the most efficient use of the 
input data. 

(2) Delete the first line of all the other files in the catenation chain. 

(3) Concatenate the files and save using a filename of choice. 

The input data is provided for entire X-Z grid range along a given Y-line which depends on 
the value of the YEXTRP parameter and Strouhal number in INPUT file WAVC. Increasing 
YEXTRP changes the extrapolation plane for the far field data and may produce more acurate results. 
The data is provided only after a certain cumulative time period in order to allow acoustic flow to 
first develop. This time may be adjusted by changing the value of FFTSTRT in INCLUDE file 
WAVC. 


The output lists tabulated values of the pressure spectrum vs. angle from the jet axis at various 
frequencies (Strouhal numbers) at both constant radius and constant sideline distance. The results 
may be plotted using any available plotting software. 

F. HOW TO SI IBMIT A JOB UNDER UNIX 

The code as presently constructed runs under UNICOS 5.1 or 4.0 on the Voyager Cray-2 
machine located at NASA, Langley Research Center. The job is submitted as follows. 
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(1) Copy all files to directories of choice. 

(2) Edit file "Makewavw" to reflect directory or compiler option changes. 

(3) Edit shell script "Makeww" to reflect directory changes. 

(4) Submit job by invoking shell script "Makeww". It is best to submit job as a background 
task and to redirect output so as to avoid screen clutter. For example: Makeww try >& trash & will 
run the job try in the background and redirect terminal outputs to file "trash". 

(5) It is best to submit multi-time step runs as a NQS queue job in a queue with at least 40 
Mword of memory available. 

(6) Compile and load files read3.f and extrpfft.f to generate WAVE post-processor 
executables READ3 and EXTRPFFT, respectively. 


G. RESTRICTIONS 

Since the basic NEKTON code, including the preprocessor and postprocessor PRENEK and 
POSTNEK, existed prior to this contract work and were not developed with government support, the 
use of these codes is restricted to use by NASA personnel only. The codes shall not be made available 
to NASA contractors. A 1/4 in. tape containing these codes was delivered to the NASA Langley 
Research Center and was marked as restricted for use by NASA personnel only. 

Code enhancements developed under NASA support shall be unrestricted. The specific 
unrestricted code enhancements are wavenek.f and include files WAVC, INPUTW, and SHIFT. 


67 


TP-499 


REFERENCES 


1. Berman, C.H. and Ramos, J.I., "Computation of Turbulence Noise," in Cpnipytatipnai 
Acoustics . Vol. 2, D. Lee, A. Cakmak, and R. Vichevetsky, Eds. (Elsevier, NY, 1990) p. 241. 

2. Yakhot, V. and Orszag, S.A., "Renormalization Group Analysis of Turbulence. I. Basic 
Theory," J. Sci. Comput. 1, 3 (1986). 

3. Karniadakis, G., Yakhot, A., Rakib, S., Orszag, S., and Yakhot, V., "Spectral Element-RNG 
Simulations of Turbulent Flows in Complex Geometries," Seventh Symposium on Turbulent 
Shear Flows, Stanford University, August 21-23, 1989, p. 7.2.1. 

4. Patankar, S.V. and Spalding, D.B., "A Calculation Procedure for Heat, Mass, and Momentum 
in Three-Dimensional Parabolic Flows," Int. J. Heat Mass Trans. J_5, 1787 (1972). 

5. Lighthill, M.J., "On Sound Generated Aerodynamically. I. General Theory," Proc. Roy. Soc. 
21 1 A . 564 (1952). 

6. Phillips, O.M., "On the Generation of Sound by Supersonic Turbulent Shear Layers," J. Fluid 
Mech. 9, 1 (1960). 

7. Lilley, G.M., "Theory of Turbulence Generated Jet Noise. Noise Radiation from Upstream 
Sources and Combustion Noise," Lockheed-Georgia Co. Tech. Report AF APL-TR-72, Vol. 
IV, 1972. 

8. Miles, J.W., "On the Reflection of Sound at an Interface of Relative Motion," J. Acous. Soc. 
Am. 29, 226 (1957). 

9. Roache, P.J., Computational Fluid Dynamics (Hermosa Publishers, Albuquerque, 1982) p. 74. 

10. Newland, D.E., Random Vibrations and Spectral Analysis (Longman, NY, 1984). 

11. Jackson, J.D., Classical Electrodynamics (Wiley, New York, 1962). 

12. Veronesi, W.A. and Maynard, J.D., "Nearfield Acoustic Holography (NAH). II. Holographic 
Reconstruction Algorithms and Computer Implementation," J. Acoust. Soc. Am. §1, 1307 
(1987). 


68 



TP-499 

13. Bendat, J.S. and Piersol, A.G., " Random Data: Analysis and Measurement Procedures (Wiley, 
New York, 1971). 

14. Kida, S. and Orszag, S.A., "Energy and Spectral Dynamics in Forced Compressible 
Turbulence," J. Sci. Comp. 5, 2 (1990). 

15. Passot, T. and Pouquet, A., "Numerical Simulation of Compressible Homogeneous Flows in 
the Turbulent Regime," J. Fluid Mech. 181 . 441 (1987). 

16. Sarkar, S., Erlebacher, G., Hussaini, M.Y. and Kreiss, H.O. "The Analysis and Modeling of 
Dilatational Terms in Compressible Turbulence," NASA CR-181959, ICASE Report No. 
89-79, 1989. 

17. Zeman, O., "Dilatation Dissipation: The Concept and Application in Modelling Compressible 
Mixing Layers," Phys. Fluids A2 . 178 (1990). 

18. Erlebacher, G., Hussaini, M.Y., Kreiss, H.O., and Sarkar, S. "The Analysis and Simulation 
of Compressible Turbulence," NASA CR-181997, ICASE Report No. 90-15, 1990. 

19. Kida, S. and Murakami, Y., "Kolmogorov Similarity in Freely Decaying Turbulence," Phys. 
Fluids. 2Q, 2030 (1987). 


69 


NASA 

Matcrval Aercnauics and 
Soace Aomnsuaicn 


Report Documentation Page 


2. Government Accession No. 


3. Recipient's Catalog No. 


NASA C R-187616 

4. Title and Subtitle 

DIRECT COMPUTATION OF TURBULENCE AND NOISE 


5. Report Date 

September 1991 

6. Performing Organization Code 


7. Authoris) 

C. Berman, G. Gordon, 

G. Karniadakis, P. Batcho, E. Jackson, and 
S. Orszag 

9. Performing Organization Name and Address 

AeroChera Research Laboratories, Inc. 

P.0. Box 12 
Princeton, NJ 08542 

12. Sponsoring Agency Name and Address 


NASA Langley Research Center 
Hampton, VA 23665 


8. Performing Organization Report No. 

TP-499 

10. Work Unit No. 

324-02-00 

11. Contract or Grant Nr>. 

NAS 1-18849 

13 Type of Report and Period Covered 

29 March 1989 - 29 June 1991 

Contractor Report . — 

14. Sponsoring Agency Code 


15. Supplementary Notes 

Langley Technical Monitor: J.C. Hardin 

SBIR Phase II Final Report 

C. Berman and G. Gordon: AeroChem Research Laboratories, Inc., Princeton, NJ 

G. Karniadakis, P. Batcho, E. Jackson, and S. Orszag: Vector Technology, Inc., 

Princeton j — Nd 

16. Abstract 


Jet exhaust turbulence noise is computed using a time-dependent 
solution of the three-dimensional Navier-Stokes equations to 
supply the source terms for an acoustic computation based on the 
Phillips convected wave equation. An extrapolation procedure is 
then used to determine the far field noise spectrum in terms of 
the near field sound. This will lay the groundwork for studies 
of more complex flows typical of noise suppression nozzles. 


, 17. Key Words (Suggested by Author(s)) 

Jet Noise, Acoustics, Computational 
Fluid Dynamics, Turbulence 


18. Distribution Statement 

Unclassified-Unlimited 
Subject Category 71 


19. Security Classif. (of this report! 
UNCLASSIFIED 


20. Security Classif. (of this page) 
UNCLASSIFIED 


i 21. No. of pages 


22. Price 


NASA FORM 1626 OCT 86 



